{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Making 2-D Laplacian inversion matrices in to order to visualize it"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import itertools\n",
    "from sympy import init_printing\n",
    "from sympy import symbols\n",
    "from sympy import Matrix, diag\n",
    "from IPython.display import display\n",
    "\n",
    "init_printing()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Creating dimensions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Note, these can be max 9 due to the current index convention\n",
    "nx = 3 # Not including ghost points\n",
    "nz = 3\n",
    "\n",
    "startXIndex = 10 # The x indices are the slowest growing indices\n",
    "startZIndex = 1  # The z indices are the fastest growing indices"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Creating the points in the mesh"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAO4AAABLCAMAAABA4UtNAAAAPFBMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAo1xBWAAAAE3RSTlMA\nMquZdlQQQOkwRCKJ3Wbvu81sHhybgAAAAAlwSFlzAAAOxAAADsQBlSsOGwAABtFJREFUeAHtnAuP\n4jgQhB3yuL0hwMzl///X67Zdbr/ipiX2tOKItHIIVW5XHJzwDVo3HX67uLfebiGlc9MxL7Stb53W\n3Tnj5eC403snlXT3Nu42L/NVFOqeVe+sBqt+UKAT97Kvh+XCtuqd1WDVDwq0cbfja/tSp1QEVr2z\nGqz6UYE27kqXt2Wz6p3VYNWPCrRx95slrHNWvdnwygJ13PvlcbvMzwe26p3VYNWPC9RxnTsWDrvs\ny84tNZeNd063oHfTt1fo+lhgn79n7lg3hALQoT0dDxKQwA+pMDRxr/4mfKesKz1mzXf6F4Kc9B/0\n03zxnwFd74Jhp1vd5ef5AugY7clo+HAoQDvrQeezNDRxw8Jw4wcPUs8UmyXnGxaSNcRV9XEhuVHc\n6bg+XQADQXs+oLRSbbuPW4yoiesXho0G4twj3H1nmoTzDQtJiMu6sT6uVA+6112PcMMbG1BAOh7r\n01K4bBy3HFET95u/KUxeeKMr2bntMXzE9HqSpbiK3sEQL7anC6SBPFlgWhE3MzRxH7xSxbi0uy0/\nfuni4N3N6+mdGFfVOxicu3HPqiHqoUPbHYw/GA07PW7w7BaGOi49wpCEP1V0MYegw6Uq6EmcZldZ\n2pLB7bjfPVeALua4ZqKlqp0tFrhfY1ySiKGO+xXOiF+e40drGj1BBz11KXGHepcMd6Slkzt4RE96\nnoSgQ9vJSodiAr4y8dkVQx13D2fwxnN8bP6qvx7+M9zvPOpjXF3vYPiiFfNr0g1BDx3a/mD80WBY\nZ9qO+V4Zirh0p/0J0RY69Xzf5fvF/REXuKaG6En94LfHeu4yFph+1nX93hSD6NEx2mYs4YAY6PWV\nZ7c0FHFv3xsemPd554eeiZ9J+HPc3UR/nW/Hha6fsd6J4eFRimYQPTpG2x0PhZME7v59XNZqREXc\nu+17PfEQEwegC8VosOrVAkXck1P2Roc/cd9oMpson9ltTskbHfjM7htNZhOlM7tWim3VD6h3Mzx/\n4IUFOnGtFNuqH1DvftwXFmjj0heoD1bvn3c+Clh1rqjesRqs+tGI2tkVNlQN8+SlVZ9Y0kl/zeFX\nFqjjWim2VT+m3k1U+lZh5PxjQx0XUDph9Ui/OwOJhyqsruprrK4aKqyu6mMBGl7A6rmhiVtiddDv\n87QlVtf1FVbXDaEA6Liur7B6aWjihoUBWJ1ZQOB0Z4GxkARWpesrrK4bQgHgdF2fVqqA1UsDx/3r\n198Sxi8MCavn9Fs0+R4WkhBX11dYXTegQMDquj4thQGrl4Z/fvH0ZeC8wep0ts5YFccGJRcSOdYn\nA9/CQsdo87Mo+yiQsHryiabYiwbB6pmhuZg9XBasTh15+l10mL+IMDoDr2M96HXW8dgQC2R0fKxH\nAWD1rFD4e1c+uwFK51g90e88ZNpPlDzN7ljPvySIZgjRpj6LHdGDjo/1KJBh9cxQz26E0oLVhX4X\no8CLRL0RV9G3WF0xpAL8oWOsruhbrJ4b6riRekesTuTb02+Ea1tQclzMmr7C6s8WEDr+ZAFg9apA\nEVegNLA66Hebk4+IHlh9rM8MEKLVCvD9hPn+WJ8VoA49Vi8NRdwMSkesDvrdH43ogdXH+ox6Q4hW\nKwCcPtZnBejseKxeGoq4Vopt1avUuw798gJF3Lra+73+xH2/OZVEn9mVc/F+e5/Zfb85lUSd2bVS\nbKv+g9Xl9Pf2Pli9d1bSMfrOePaHgfZiBnxKbmXHqie2oPRYvW3Vjwq0cYUNVWVPXlr1iSWd9Ncc\nfmWBOq6VYlv1Y+rdRKXvNf8pVsev0DsDiYdKrK7rK6yuG0qsrusrrF4Y6tmNUBq/Vsev0M/Tllhd\n11dYXTeUWF3XV1i9NDRxw8IArE6UAr+jO0mMhQQ6tCdyxq38ltBuzRD0wOr6gNJKFbB6aWji+oUh\nYfVS3IuAhQSjRtvT+mPBILRbM6AAfq2u6dNSiF+r54Ymbo3Vc3EvAqg3dGh7Wn8MBp5nxuqaIek9\nntP14PYJq+cFmrg1Vs/FvQg1Vtf0oN7UV8DjmqHG6poeBRJWzw113ECxM6yeiztxE/WGDm1H6w8l\nA36trhhEH7G6om+xem6o4zZYPRd3IiTqDR3ajtYfSgbQbsWQ9MDqir7F6rmhjltjde2jVWN1Td9g\ndc0QCghW1/SxgGD1wlDEFUwOrE4rg/8Ven+uRC+6kT6j3kK7RwYpAKwuhbpDEgO97bF6aSjiCian\nD5b/tTpwebfvDGJDh7avzwyg3WODDAhYfazPCgCrl4YirpViW/UfrH5yGfyuw8Xs/q4if06/Ie7/\n6r8J2fg/0FgW/FX9z5mLl47E/zchy+L+BZ1IYMDY0U1uAAAAAElFTkSuQmCC\n",
      "text/latex": [
       "$$\\left[\\begin{matrix}f_{03} & f_{13} & f_{23} & f_{33} & f_{43}\\\\f_{02} & f_{12} & f_{22} & f_{32} & f_{42}\\\\f_{01} & f_{11} & f_{21} & f_{31} & f_{41}\\end{matrix}\\right]$$"
      ],
      "text/plain": [
       "⎡f₀₃  f₁₃  f₂₃  f₃₃  f₄₃⎤\n",
       "⎢                       ⎥\n",
       "⎢f₀₂  f₁₂  f₂₂  f₃₂  f₄₂⎥\n",
       "⎢                       ⎥\n",
       "⎣f₀₁  f₁₁  f₂₁  f₃₁  f₄₁⎦"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# The function at the current point\n",
    "f = []\n",
    "for z in range(nz):\n",
    "    f.append([])\n",
    "    xStart = startZIndex\n",
    "    xEnd   = startXIndex*(nx+1)  # +1 due to ghost point\n",
    "    # +startXIndex in the range as the range does not include endpoint\n",
    "    for xInd in range(xStart, xEnd+startXIndex, 10):\n",
    "        ind = str(xInd+z)\n",
    "        if (xInd+z) < startXIndex:\n",
    "            ind = '0'+str(xInd+z)\n",
    "        f[z].append(symbols('f_' + ind))\n",
    "\n",
    "mesh = Matrix(f[::-1])\n",
    "display(mesh)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Do the vectors"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAC0AAAF2CAMAAADN+uEoAAAAP1BMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAADFBd4eAAAAFHRS\nTlMAMquZdlQQQO0wRCKJZt3Nu+98bODTYm0AAAAJcEhZcwAADsQAAA7EAZUrDhsAAAZmSURBVHgB\n7ZzpmuMoDEVJvHRPxUklPX7/Zx0Jc0HCgDXZvkyN+VF4ucFYwImu3Wl3mH05ula5LCLnDnPXUxla\nYndmyXFm9aEpTCfP/0P1SEE6ffVTn8KgtlRMxsm586V3l2/S0Gdol8rh6iv+o9TT6Nw0j2440QmS\nDjS8h+54KatJ7L5DSxceMfokfaai5jbmpc/jTO27G8+GhnqZBSPNBm72cm6oh9H1NBNcH9V8oUrb\nX/PgrqQ+UYuHpScN9Xjs+7Gfznz50c/L+aveNp1J5cLCzZhA33dLvOnvDcf06MSjvDF1U0dxOXWX\n+YiZoMZSqUs7u3odlT0mWzF5nD5Td+WZ64uKd4k+E5HiyOziotTcRE6fC6mXJR3Uv+Zf/MGl5PS5\n0fo8+dVM5/9k3yQF+tDK9KuZ1KonZfoQskrrskIfWs+hl6rtCn3OUazUaMHXkT5fBPKv8IWq+i3l\noM/hexiGaxieqhr0ufnv99BOXS0vhO3PUo/0zWErJxrUPT/RsXrfWD6BPv10xFirfpfo09G3fodc\nSalL9Ok41+I0hAtvNOnDog5ks9BnvCErVj0p0mfsvwEfRYgafcp3WaEPL4CQyaue8B2hBPr4Pp9m\nTp6oVNWgD/P7fLPS50Dp7JEzSS7VtpfT2d/PUu/00cPzWaOTaH8/faLjsjgvOC7UPjgqJpI+yedk\nfqdEH2oKKtR0qEafslr1RNBnWy3ps62W9NlWk0IW3B1qOqf6LbXJcVmcFxwXat9So219pQ9Uq3lS\n6G06lM+TdKa09V+OSVrzpTuTx953l4/Tx+S8QJ2V84rfUiX6KOe1ypGD86K4LetROa88gnBeUc3h\n3nReSr3pvKR623kJtcV54S6dyXmBPjbnBer8GOe104dmlyj5ahCnCpvvUz+BPhbnBfqsnFeTPsp5\nbdKH4xidVx7BEn22nRe1GDINi/NKau5JyV9Wch+L8wo9MTov0MfmvECfH+O8dvrQDBQlXw3iVGHz\nfer76RNzHovzQs4DCvl75rss0iflPMJlVOmTch6hziMo6EMrs/nuSDsvPG2utJ05Lzxtrqgz+iDn\nqaj1xIg5j0Wdch6L80LOAwrF0fnEbF1HZb3HI7/TR8clXw36bL73PvXj9IkUoptQ/ZZPnUEf1P6G\nWd2kT6IQ6TfpkyiU94T2C/QBhXJ1mT6gUKau0AcUytQV+oBCmZp2RQF9UPMpFW+hdaAPan+uqgZ9\nULfV8jpxu9p2VMiNV6t3+shoN+aJloW9V49O+k57An0szgvUsTkvUMfmvCR1LM4rUsfkvAJ1jM4r\nUcfivBJ1tt95hec9RucF6ticF6izO6/iakwH37cu0zXLW5/ak/t5UrpPdZcym3H8dta/yq05L/kc\nGdRZOS/1LYXnyIk6wmXkbydUNrNQR6hVv7NsJlCnotbZTKRORZ1lM3jeU1GvIrxQx6AW1DGoXaKO\nxXmBOh/svGKOvAphfmCVI+cCta/miTpT2tnV66jomNxPH/DE5KXAE2Q1vluqJ4o+4V/uIKuJ6iJ9\n+CzzRGY1rXfoLmYxhrdYkSdGL4Xnx6a3WHh+nLKaqjtKPDF5KfAEWU2MYMpj/aHlD3iCrKatFh9M\nm2os0+HK1qvVO3104F8d7+Ks0l0Ie7onT6APshpqX7Utcx/QB7Xvi1KX6JOymtB2kz7cZvRSrdwH\n9EGd9VvlPqAPar6I7LfOfbgH4e0VaqVe5T7wUKiVmq+FAvqgXo7LnkDpa9AHtT9YVYM+qNtqdSHs\nVNuGQNWvVu/0UeGuziqtwt6rRyex6n76UF+XhVbzUpI+JPaJxspLxZ5I+tDjZf7JF39I/3KrQp9+\nXKnr9DkMa7WKt6LPRL+wKvUE/Vb0OZ821JI+I/1r6HbbHACUoaMyM7zzmKAnUKI+hX5bfrnlztf5\nOPycX27t9MEsWGo1v/Wpwt771I/Tx+S86B49fWzOC/RZO6/ff6Vo4bkP/8KTV5pyXn//xpur5QPL\nGmWORPqELvF5Fe8yfSzOK9EnPU9WbVfoY3BeiT4m50X35Oljc14UAE+fH+O8akTmsdZFjbw+Vdjb\n1eugPBKTJ9DH4ryo054+NucF+qycl5pVOX04MNF55TEp0cfivJD7mJyXoE90YBb64PkP3UDeb74p\nLqCP0XnRJzx9bM4L9Plo5/Vv/l8LmtVc+If+9eL/X4u+d/8AExCdxZ4JfPkAAAAASUVORK5CYII=\n",
      "text/latex": [
       "$$\\left[\\begin{matrix}x_{0 1}\\\\x_{0 2}\\\\x_{0 3}\\\\x_{1 1}\\\\x_{1 2}\\\\x_{1 3}\\\\x_{2 1}\\\\x_{2 2}\\\\x_{2 3}\\\\x_{3 1}\\\\x_{3 2}\\\\x_{3 3}\\\\x_{4 1}\\\\x_{4 2}\\\\x_{4 3}\\end{matrix}\\right]$$"
      ],
      "text/plain": [
       "⎡x₀ ₁⎤\n",
       "⎢    ⎥\n",
       "⎢x₀ ₂⎥\n",
       "⎢    ⎥\n",
       "⎢x₀ ₃⎥\n",
       "⎢    ⎥\n",
       "⎢x₁ ₁⎥\n",
       "⎢    ⎥\n",
       "⎢x₁ ₂⎥\n",
       "⎢    ⎥\n",
       "⎢x₁ ₃⎥\n",
       "⎢    ⎥\n",
       "⎢x₂ ₁⎥\n",
       "⎢    ⎥\n",
       "⎢x₂ ₂⎥\n",
       "⎢    ⎥\n",
       "⎢x₂ ₃⎥\n",
       "⎢    ⎥\n",
       "⎢x₃ ₁⎥\n",
       "⎢    ⎥\n",
       "⎢x₃ ₂⎥\n",
       "⎢    ⎥\n",
       "⎢x₃ ₃⎥\n",
       "⎢    ⎥\n",
       "⎢x₄ ₁⎥\n",
       "⎢    ⎥\n",
       "⎢x₄ ₂⎥\n",
       "⎢    ⎥\n",
       "⎣x₄ ₃⎦"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAACoAAAF2CAMAAAAvJvpRAAAAP1BMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAADFBd4eAAAAFHRS\nTlMAMquZdlQQQO0wRLvv3c2JImZ8bDIaVl0AAAAJcEhZcwAADsQAAA7EAZUrDhsAAAZTSURBVGgF\n7ZvrsqsoEIUxGs+ZnYvJjO//rNMNrr4IiJWqVCW79EfwsiQtwpdeREM3x+UUqsstKULo5n6g5VxV\nhicfP80s7eoqPfI00vFyvemRbM1KQ7hvxBu8dN6I10t/5jH7Wt3hap22QvW13k/D0FdjsLWO850C\nmB8hDNMwxW/uLuUAzjHUeQpP0p2pMbr+ZEKytfZ3qmIk6Y1vSTzvXJHee1L8zOeRYwhXDromnZ90\nsJ/HLlZ4462q9Ie+//qkPsHNexs2pCeq50JBdCmALenYD6ymK4uXRd9RDYAOpeXGqs0WgDIMFAW3\nK31eZae/sbJ76qeeLu3R3+YThxwXewuwr1Ie0hdb4EBGjoypv3BnTItt1wwZEw3yE8MhLlaaIeNG\n0jQmWcvSP/OfeFqGjCuNscfMA42Xfw3gC8ig0RWHI0ttAHz+ChkLOVjppAVkhIk5lhZbawEZT1W6\nWnF6EGT8EGh/8KNmaxUpkNHdz+fzBfegKKUIIzKu8acVNZSlOOrKj5COAy7SxbbeeBATX8gH1tWs\ntl9sgQMZBWQM00nupG3XDBk9/eL2kpJYaYaMnsYLC9LCa1VksCSeH7UtZIxXjEI3YnNkjMNdfo2d\ntISMymUVkMGdGemcbYF0ofS5ICPG+YjQ40NFKZDBfH1ecQ+KUiCjo9TwxOlZXMpSHHXlR0gPZCgJ\npRe527RsuLv1OcYkwJCg5HBtrIIMGBKU6bqsVJChmfvaQqyRQZVAgpLrLSGjIrUBABk7pEDGDimQ\nsUNKEiy4HJS838YKHZUwJCirUhgSlKmOSq3mC2T1I6ToAxJVecX2gbJC9r7/so4sQ7MMIKNtTICK\n3JjgZyNDhjcmNtOEMaG7nsaUNya2D8CYiJR7yrYxcdI4pcHnuBFbQEbLmJham8YElxXaxgTI2GFM\ngIrvNiYHMn5HlsGpLpJiN7ayLCM3JlVkeGOyiQwe0mpMWshoGBMd3E1jolIOoOi3CllG05gsyNhj\nTICMHcYEyPhuY3Ig47uRIdlF05ggu0C2QUOMFh7cGTI0u7BpeREZml1YaQUZNLo2/4UwxgTTnrVa\nTZaB7KImNchAdlGTpjbhT8kumlLNLprGBNkF0JG+zjaWBlBc+wjpgYzfgQxBB3c127MkywAyUKY+\nydIqMhQdLN5EhqKDpTYAY0yADC1X0uUfE64i/cdqSi8tIAPo4JNtAAVkAB1rKW+nBchAGffaWiEM\nQAbKdKAoBTJQbkilerdSrNUpZON90gMZvwQZTWMCVOwwJkDFDmNiUdE0JoKKtjFZkLHHmCgqmsZE\nUdH4x2TJMvYYE6BihzEBKg5jImSLKy/i7ZjTNHOavkXTlm1XSR2CzGHUjInMaQIVuTHBz4bMaSoq\nbFpup8Bd6pBQYaU2Vus2FlTUpJo6CCpqUpM6YA6jJnWNmVDRkhpUtKRBUdE0JkDFBxoTZJqusfIN\nm2nmR90e2wfcgXzjRemBDEUGONB2G+AAUoh0O+wtUGQsD1kghVBphgw+xOfZFKL6zyk9fLU8ZJFA\nxifbABQZwgF1HV6qyOAA0lxm4/ks/rblIQtNIXytUUIfyoG22wAHkEKkOuxlodYADiCF2JDKOW6l\nWKtTyMb7pAcyvtxtLC4D6IhdxnYXyTKADJSpc1lphgxNIVjM0ioyWJBSCF6rZRlABkqS2gByZBh0\neGkJGUDHqlaXZcBloPS1cuy8ABko014X67KLkLY8l4UyHbCXJVIgA+WGVM5xK8VanUI23ic9kPHd\nyKAukn5ia25DkEHKuJ67DbSAIIPGzJQe017n2gVkDGMuLSOjo/dM4r+761oRgCJjotc8NqWCjOej\nJQUyRno4tFErBty5p2VmuMrjZ3yoNgwfS6zNF0fC8zKfzl/+4siBDPRXercIDwaj59jSdZcPeqRz\nQUbbmNDFRGTsMCZARsGY/P0nNorMZdDbbnEYemPy3199aUDnMgQZMrdJddl2LSCjaUwUGTq36Wst\nIKNlTBQZbWNC1xCRscOYUFwRGd9tTHTAxK5Q+7A3tqZZ9r9fesxl6FwGNXrMMtrGhJQRGTuMCZCR\nGxN0lwwZfPPVmNg+UEJG05ggy2gbE4OMhjExWQbmNDhqGytv8wJk7DEmJI/I2GFMgIyPNCa73z0f\n+c3yYcB7hbHB/Ed893wYwv9BCawgcASXkAAAAABJRU5ErkJggg==\n",
      "text/latex": [
       "$$\\left[\\begin{matrix}b_{0 1}\\\\b_{0 2}\\\\b_{0 3}\\\\b_{1 1}\\\\b_{1 2}\\\\b_{1 3}\\\\b_{2 1}\\\\b_{2 2}\\\\b_{2 3}\\\\b_{3 1}\\\\b_{3 2}\\\\b_{3 3}\\\\b_{4 1}\\\\b_{4 2}\\\\b_{4 3}\\end{matrix}\\right]$$"
      ],
      "text/plain": [
       "⎡b₀ ₁⎤\n",
       "⎢    ⎥\n",
       "⎢b₀ ₂⎥\n",
       "⎢    ⎥\n",
       "⎢b₀ ₃⎥\n",
       "⎢    ⎥\n",
       "⎢b₁ ₁⎥\n",
       "⎢    ⎥\n",
       "⎢b₁ ₂⎥\n",
       "⎢    ⎥\n",
       "⎢b₁ ₃⎥\n",
       "⎢    ⎥\n",
       "⎢b₂ ₁⎥\n",
       "⎢    ⎥\n",
       "⎢b₂ ₂⎥\n",
       "⎢    ⎥\n",
       "⎢b₂ ₃⎥\n",
       "⎢    ⎥\n",
       "⎢b₃ ₁⎥\n",
       "⎢    ⎥\n",
       "⎢b₃ ₂⎥\n",
       "⎢    ⎥\n",
       "⎢b₃ ₃⎥\n",
       "⎢    ⎥\n",
       "⎢b₄ ₁⎥\n",
       "⎢    ⎥\n",
       "⎢b₄ ₂⎥\n",
       "⎢    ⎥\n",
       "⎣b₄ ₃⎦"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "xVec=[]\n",
    "bVec=[]\n",
    "# Do the inner loop, so start ranges at 1\n",
    "# (nx+1) to include outer ghost point, +1 in the range as the range does not include endpoint\n",
    "for x in range(1, (nx+1)+1):\n",
    "    for z in range(1,nz+1):\n",
    "        xVec.append(symbols('x_'+str(x)+'_'+str(z)))\n",
    "        bVec.append(symbols('b_'+str(x)+'_'+str(z)))\n",
    "\n",
    "# Do the inner ghost points\n",
    "# Must count backwards since we are inserting in the front\n",
    "for ind in range(nz,0,-1):\n",
    "    xVec.insert(0, symbols('x_0_'+str(ind)))\n",
    "    bVec.insert(0, symbols('b_0_'+str(ind)))\n",
    "    \n",
    "display(Matrix(xVec))\n",
    "display(Matrix(bVec))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Make the global index matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAxoAAAF3CAMAAAAVa83NAAAANlBMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAABHL6OuAAAAEXRSTlMAMquZdlQQ\nQN0iRInNZu+7fNewZkkAAAAJcEhZcwAADsQAAA7EAZUrDhsAACAASURBVHgB7Z3rdhwtzGbLh3iS\ntA/T93+zw1GPAAGP3lnzrVkp+odNVbERqITd7nTvXC/P9Hi9zuNk4GQgZOB33hHX9fJ8ew+PXycr\nJwMnAzEDX3E/vD7j1ng5GTkZOBloM/B1tkabkHN0MpAzsN8aH2+P9zf/r5XPD2eK/7y9/fw44/x5\nPB4//+GJ4Mubb3Kf7x/Xy9tfH3Q9Hu8Pz4Jefr++hQU9Hp44H+/vD+/M0i3944lS7qarFlABaO1i\nlp6uWqiMpxbUjKa1sN8an6H0Pn67Enldf36ezq3xEV8FeDxd5fcnVviv59cu4cP13z/DqeWJ7/hn\n2fuyy3Dx4zOs5e1zOD8/8ZX/+Hu6Xg9Jm/zV9ePh4ztk7M83f0vlbjpqQRhHLVTGUwuV8dRCZdKt\nmNbCdmt8fccBHp5bfH38vL16t8Zr2krfKdq8etor789YEs/f7dn90fu3c2t8Pl7f+UrKE/iMP/x/\nPIHecgjXzviVfi78cd2e15TkV3ZquJt8LYBBa3df0JOvBTB8LYCJM5rXwnZr/KQ79ctb6e9e4DsV\n+KvrFYG/3+EH88fTVRchG3++PtmyKLfTVa2Z+Xp699KVn0i5noNdj7SSD9dPh5ywh+MOlbvpqgVU\nAFolndNvpaerFgrjqgXMaFEL262Rf73/TT+fp0saLyD4eM0885nulG9rpIH8T6ge1//A1vhx/f5D\nSv76nrf9ff6EX7cPz3PKj2faTu+OJ6/lbrpqARWAFpZpt0pPVy3o0dlaALOohd3W+HimJ7N/vc+0\nEdxOwuTsb8ePsjLEp/NP6uvrj39rfH15X4p4fr6Ely9cfzqlBTl/n12vz+eXa2dc13f6Nftw/I2W\n76avFlABaE3uupxuepK1oBm2FoRZ1cJua/x5pl/zL/mbrGHbkODbnrrDX2+Yl8dv7874CD+Wvb81\n4k/mP46fsulpXvxJ/u36+zgAzioPxNvzqV5v0dmctR/hBfuQA8dPu3w3fbWACkBrNqN6Xvdka0EY\nRy1UZlkL/39tjemrBTV3xvfHt+cF0lB8YQzv1khhfzxP6D+e6fff27fvhbqPVLXGKqenvt5ePp/x\nby7H4zP8dPjl/63xP7s12FqoZR7Xz9ZCZZa1sNsaH/nn+P/ME6o377OJVA7frvL7FTfSf9oab64/\nrPPrZl+uXzXh9RLvSwpf8RWCr6crBemNEC8Px5+PuZR8tVDLLyyKfpasetK1oJjwW5pLRGHWtbDb\nGlf+W+OXI4+pXpsJpzP7L+/eF4L+pBeBfjxz+0jPD51b47U8Off8bM5P6L8cT+hjgtgflZLMXAov\nnhQU1vPyermbrlpABaAl05400JOvhcK4aiEzm1rYbo18s77ojV8WjUVOsjCe/hX/asgrHC+aZ/KT\nlh9P+f16jY/n96vnhaBc5m+uHORXqJy/NT58/9wX/6TJiXn1vESVkU/HD6JyN121gApAy7yL6qT0\ndNRCYVy1kJlNLWy3xnt6EdL1r7pxrbJItfB1M/+Dvetf1vK/3P12PdFJk3D+k1/+S/+368lO/peg\nN99LuOE90OskDVfLv2n/eP4R5StCL56slbvpqgVUAFrD9LsTtaenFgrjqoUaJ4af1sJ2a1y/4xtF\nnH/q/oet8ed3eA/M26vnT93rK/6s/OV4paXeifzKfj3afk9vQkgFte2KDj/hqduH8xWqX+6t8Te9\nOvXl+R0Y/qD5uD4+Pb9nail5aqEynloojKsWCuOqBcwtvJli9gfufmt8PMJb130vAl1vr9/P375f\n8p/5DUSurXH9fQvvSfT8BZAq9/HzfPom9+ft8frme60pvFzy9ubN2x/nVgrLeQlvSfT+80mYmSdr\nuJt8LYBBCz827Jb0dNSCMI5aABNu0rwW9lvDXsY5ezLwj2fgbI1//Aaf5f3XDJyt8V8zd7h/PANn\na/zjN/gs779m4GyN/5q5w/3jGThb4x+/wWd5/zUDZ2v818wd7h/PwNka//gNPsv7rxk4W+O/Zu5w\n/3gGztb4x2/wWd5/zcDZGv81c4f7xzOw2xrls5UuR5l8HtPhKMuMT1JW4rgkZZWJ3jnq3aqwhdF+\nMiCheCQZ60ICg9aaCG/gF6kdr6oDE0ef6sma0GD4MgAT3qdEquqEcZSBMBdfBopZlcF6a4jLyuEo\nE8bhKKuMR1JWmSu9j5uTlFWGl5QpWxjrJ1MIrScDg1ZTnsYBevJ6MjBpQOqTU4qhy0AxdBmA4csA\nDF8GYNZlsNoaymVFO8oaJqSfcZSB4SVlYHhJGRheUgZbGO0nA4KARmU3p8Cg1XQwDtCT15OBiQPO\n9WQ6nGLoMmiYMBhTBmD4MgDDlwGYdRmstkZMXvlkm+MTYZVxOcpKnPTp1IuUlFUmvd+elJQVhpeU\nwRZG+8mAqATqWjPaYNAyujWn0JPXk4EJQy30ZDqQYugyAMOXARi+DBRDlwGYdRn8v9saLkdZ3YLx\nlrCSssK4JGWZcUjKYAuj/WRAwmL0wnS9dW0waHVdhsOuJ6Una5iFnkwHUwy9NcDwZQAmRafKAAxf\nBsJsyoDdGg5HWakGl6NMV9DsU1f6dsV2ZTySssJ4JWXBFubzk8X/2Cd99KlOsp+8eQwpGVpmR3VS\nerJ6ssBmZqUnUwFKMzGvjjIocVxlIHMLDbYMKuMpg8Ksy4DcGh5HWa6GoKJ1OMpUBdGSMmEckrLC\nOCVl0RbmkzBdVTAmkwx3Y/eoTPjFSavqSk+HnqyOvtSTDVPNcTxlkOP4yqDOLYSny0AYRxkUZl0G\n5NZIqSIdZbkafI4yVFB1ZAw3ZzhRGY+krDI+SVl8Ice5NeprPzXgMHvjRGU8zh0wrJ6sjh6f0PPS\nIcS5yDLIcXxlUOcGVYqRpuFUnpunDGqcZRnErfHxa/qx5+bOko6ywrgcZYjDS8oK45KUSZzwv7XR\nkrJkC/P5yUQwJgGHGzqcEOZCa+jUndA9ST1ZGX2tJ+vCNDMiy6AwrjJAHL4MCuMqA4mzKoM/wRey\n+L/8yp11OcoK43KUoYLUj6f+/nTHNU7a16SkDHHCYKSkrNjCPH4yCMaagN0C2kPF0H/tFsanJ0uj\nb/Rk7dSuEsdXBnkVvjKQlTvKoMRxlYHEWZQB94Qqr490lJVqcDnKpIIckrLM1CdgnCJE4sR7z0nK\nqi2s/NpmPG0VCTGagDHo7AEGrVnfer729OjJCrPRk9UI+XuN4ymDynjKoDLx+ZQu3nY27VFhXGWA\nOHGsWRlwW8PlKCvV4HKUSQU5JGWFcUnKCuOQlIktjPeTCRLyLgtrb+hwBAatoVN3Qno69GTCpLGm\nerImkjCOMhDGUQbCxCcypKtOGEcZCLMuA25ruBxltRo8jrLKBOEamROpOpekrMThJWXKFsb6yRQi\nk2xKzTgAg5bRrTmFnryeDEwaaqon04HA8GUA5qLLQDF0GYDhywDMugyWWwMuK95RBoZ3lCmGlpSB\n4SVlYGhJmbKFsX4yhSCgrjWjDQYto1tzSvWkVXWKWerJdCDF0GWgGLoMNMOWgWLoMlDMsgyWW0Mn\n6LRPBu6VgbM17nW/z2rpDJytQafqdLxXBs7WuNf9PqulM3C2Bp2q0/FeGThb4173+6yWzsDZGnSq\nTsd7ZSBujf/1/F/3WvRZ7cnAPgP/O2yNxdsL9wOcHicD/2YGzhOqf/O+nlX9X2fgbI3/6xSeAf7N\nDCy3BsxftJ7sAhMSxvrJwn/++hM+OKKtY7t0Iw5aPBN75jdg8gztJ9MTovVkNQUePVllPHoyMCs9\n2ZCVnCy+DJoEk2UABhK1YSLGiTw3nXWjU3cKNx+trstqayjzF6snC0a9EKDILaoQrQs5HIKBPWvo\n1J0Ag1bXZTjselKfllEM6ydTCK0nA8PrycDwejIwaz3ZkLucLLoMEl8SzJYBGL4MwGBlw9ytE7j5\naHX9VlsD5i9aTxY+nxB/+sdPPfJ+MjCwZ3XTHA7BoDV06k60PTk/mWJYP1mDhCkwejIwvJ4MDK8n\nA7PWk3WZKzI3vgwinxPMlwEYvgzAYGXx3O6Bm49Wz6y2BsxftJ7sAhMiyacw+qjtMRjYs9oe4xEY\ntMZe7ZmmJ+knUwz7uTMgvJ4MDK8nUwytJwOz1pO1iasyN74MAo8Ek2UAhi8DMFhZN3frEHNDa+i3\n2hqpc3pyROvJ8vjVFsbnpDwJE3vWME/zRI0jT+HMXu3JypB+sgRnht0aQHg9GZjUovRkYHg9mTAb\nPVnqhy8lWa4yQIL5MiiMqwwQhy8DMGhhtaW13RrR/OXVk1VbGJ+T8Bc7Pt0n1rFhtt0JMGh1XYbD\n0tPlJ8uMy0+WEKeeDIvg9WRlas/nFy9uSkz+nPfjGWVhu0dJlqsMVILpMlBMmBJXBg2DDC6XBAat\nEVhvjWL+cjmYlC2MzolixLg1zrU9Awattsd4JD0dfjJheD9ZQVx6MgkTZs1WORheT1aZtZ6sSV1N\nlqcMKhMHYstAM2wZKKaurJm7dQAGLaPfemsEIJq/PDmJMaotjM2JZqo9K57bPWocRNwRtafPT4Y4\ntJ8sIl49WQ1T7Rj7xdTl+PRkOc5ST9aErsnylEFl4kBsGWiGLYOWqRlspj8cgEFr6HRdcWt8vE8V\nbYEI5i+fniwzMRabk9g3xknftXUsnVh8qQzoRedyKTI+P5kenfWTpeX49GQSxqEny4xPT1bjrPRk\nOo2SLEcZCBMHIsugYUSipmcytlumrmzsp8+AQUtfr+0/7/P3UMH8xevJwIQIZE5ahvtrFwxadU2z\n79LT4ScT5qL9ZEB4PRmYMPnpK+3twsDkHyqMpQ5MHouw1Klk0WWgmBCHK4OOocpAMf3K2lypIzBo\nqctorp5QwfyVb9YXoScDE2JwObk009qzMM++BQatvk9/LD0dfjJhLtpPBoTXk4Hh9WTC1CdghKVO\nmJKbmZ5MpU4liy4DxYSRuDJoGLIMFNOvTK2gbYJBq+1RjlZbA+YvXk8GJgTgcnIpJv+r/fv+f9kD\ng5a5QHWy65kP1XWrCSa/gvb70+rVnAPC68nA8HoyMLyeDMxaT9YsKB4kkC+DxOdgdBkohi4DMFhZ\nOrf7UudWVmZ1X20NZf5i9WSXYuicgIE9y5qrPgcGLX3danc9KT8ZmPROhFRQ1tjqHBBeT6YYWk8G\nhteTgVnrydRycjMniy6DBNUEkz8hwfBlAAYryxPefK1zC+/cmL1SvtoaF8xfrJ4svOb2Fv5f0r9x\nYrSfDIyyZ22WBka1eCa8tPPzfBLPQNTotJ8MKaD1ZCrMxerJFEPryRSz1JN1mazJ4ssACebLQBhH\nGQijVtZN3jis6wFtdFpuDaP/OXUycJMMnK1xkxt9lunNwNka3oyd/jfJwNkaN7nRZ5neDJyt4c3Y\n6X+TDJytcZMbfZbpzcDZGt6Mnf43yUDcGkfRdpObfZbpycBRtHmydfreKAPnCdWNbvZZqicDZ2t4\nsnX63igD262R3wRJ68lS6qr1ivWTRSgxDj+ZMA4/mWJoPxlWTvvJgISApJ4MDK8nU4yI1+ISVw8w\nsVe9USsiLCF81O3lLb4vDq010fZky0BGd5SBMI4yUMyqDLZbI79Hn9WT5YxlhvaTJSgxvJ8MDO8n\nA+Pwk2HltJ8MyEXryYRx6MmEcejJhEEyUmv1BQxaq/7xGnryZSCMowyEcZSBMOsy2G2NYrBi9WQp\nY8KEI8ZPFqHM8H4yMLyfDIzDTyYr5/1kgjj0ZMI49GTCOPRkwiAZsbV8gEFrCYSL6PkZP4BNlYEw\njjIQxlEGwqzLYLM1qsGK+jhiSVdheD9Z4AoTsxjeKPySvu2+VIb2k4UBC+Pwk8nKeT+ZICEg+1kF\nYRx6MmEcejJhkIxdni8waO0g6ekoA2EcZQCGLwNh1mWw2RrVYCWj7VISrhfG5SerceLwrJ+sMC4/\nWWY8fjJZOe8nEyQsxr01HHoyHYfVk2lGJz0mfvYAg9asbz0vPR1lIEwahCsDYRxlUJlNGay3hhis\nHHqyynj8ZJVJOZl96qpmvXwX5pX3k1XG4SerK3f4ySoSJ0pvja+v9zf8tqT0ZDpOeAoDx12XKX2o\nmJoMfdlsg0HL7KhOSk9HGQiTxuHKAAxfBsKsy2C5NWCw4vVklfH4ySqTUkL6yRRD+8mEcfjJ6sod\nEqaKxOWwW0MzrJ5MMbSeDIwkI2V99QUMWqv+8Vrt6SmDyqSxyTJQDF0GwqzLYLk1OoPVT/Cfbx+V\n8fjJKhMHr3qMXSAwvJ8MDO8nS/MIK3dsjYrE7+zW0Ayt24mQ3BROT4Y4SEY6t/sicRBxh8SenjLA\n3PgyAMOXAZhlGcStMVO09QYrRk8GhveTgQmTJv1kYHg/GZigpnt/f3mk//AgJWrzJazc4SdLg5Vk\nubZGTbDHUlcZTk+W1xkZnYzN6vV6QhsRd1zsyZdBHq2MTpYBGL4MwKzLYKFoUwYrWk+mGNpPppgw\n6/xvInn686+Kof1kiskDE34yJWaj/WRNssit0TL178T5+uMVMLSeDMyQjHksxEFr3jtfQU+6DDC3\nOARXBoqhy0AxZa4TvdriCZUyWOX1vU3GyBHSV8XQfjLF0H4yMPUJ2N4OAqZMmPCTKTFbvlmEpq5J\nFrk1NEPqydTUaD0ZmCEZ6iZ2TcwNra7LcIiedBlgbmG0j/x64DBuf0Li8GXQxInjzcpgsTXyNJLL\nitaTKYb3kyWoOLNe8g/nPM7ua2J4Pxni8H4yrJz2kwEJAcmtoZj85g3CUgcm5+53ePqye4BBMnim\noxcgevJlAIbX1IHhywDMugy2WyMZrGg9Wc5Wtl79hD/0Pr7T/165SGK9VExZvzxbIzG8nyyFSgzv\nJ1MrZ/1kCqG3BhheTwaG15OBQTLqHZh9B4PWrG89r3rSZaAYWlMHhi8DMOsy2GyN6rKi9WQhN5Xh\n/WSKof1kYHg/GRjeT4aV034yILyeTBiHnkwYh54MDJJR63n2HQxas771vOoZcv2Kf7GpHYzviqHL\nAAxfBmCWZbDZGsYCzqmTgVtk4GyNW9zms0h/Bs7W8OfsELfIwNkat7jNZ5H+DJyt4c/ZIW6RgbM1\nbnGbzyL9GThbw5+zQ9wiA3FrHEXbLW71WaQvA0fR5svX6X2bDJwnVLe51WehvgycreHL1+l9mwxs\ntka1a9F6spC4yoQm6ScDw/vJFEP7ycDE+8v5ycDA7BXp5QMpQGsJhIulp0NPJoxDT6aYlZ5MT1bN\niC4DxYShuDJQDF0GmmHLQDFxmdMyWG4N2LVoPdkFhveTCePwkwmT3kj56/ml7+akLUy6zn1aBoyY\nvSajy2kgaMnFSUN6OvRkwjj0ZMKs9WR6lmpGdBkohi4DMHwZgOHLAExa5bQMlltD7Fq8niy4uULA\naORy+ckK4/KTFcblJytM+FakcLG1fMh6lHVsCYQfkjUFqrVB0NOlJwujxlS79GSFWevJ9HQxI74M\nwPBlAIYvAzB8GYCJqyw+Qb3g0l5tDdi1eD0ZmBh2/7HAOA0wvJ8MDO8nAxOCFllbjL96KIb7WKpe\njoJXMcI19IzbKjzrId7FrZgkpvkgrBZg1nqyNInyBTPiywBMGIQsAzB8GYDhywBMmNuiDFZbA3Yt\nXk8Ghs8JGN5PBibdQOoJVcOQfjLFsFsDCFppkosvXU9KTwaG15MJs9GTDVNNM+LLIPF1FeTWAMOX\nAZjUosqgYRZlsNoaYtdy6MmEiRMgc9IwAWP8ZB1D+ck0w/rJFCNmr5Tb+RcgaM175ytdT0pPphha\nTwYmf6r6Qf2NFuYYZ+Qog7SqugqyDBomHDBl0DFUGWhmVQaLrQG7Fu9gAhMnwOWkZTg/WcOQfjLN\nsH4yzYjZK+V2+gUIWtPO5ULXk9KTNQypJ1PMWk/WTzjNiC+DhMsquDJoGa4MGoYsA80sy2C5NdKf\nCm/fH3xOWiMXl5OW4TwrPcP4yTQTn29+1p9qKVn2F82kHjCV2UD82VrThtasbz3f9qx2jHrV/q4Z\nVk+mmaWerAuZZ8SXQcSxCq4MWoYrg55hykAzyzKIW2OmaBO7lkNPJkycAZmThrk4P1nLcH4yMLyf\nDExcEOUnA4JWhudfm56kngwMrycDs9aTdTPNM3KUQeCxCrIMGoYsg5bhygDMugwWirYi7PkKT0dp\nPZli4gy4V6jy094YJzLcX7tgeD+ZMA4/mTCD2StO1nwAQcvsqE42PaevtCsgNMHwejIweSzGUhd6\nlhnxZQAmtMgyaBmuDMDwZSDMpgwWT6iu/HLG1/NvSQ2hJ1MMnxPECS/RR0tQXmZozB9geD+ZMA4/\nmTCD2Ws6NyBoTTuXC7onqycTpj512Vvq2tsTYs/0ZO1864zyDmHKID6fkuKmtwYYtgwQhy8DYTZl\nsNoasGvRerIgEPoIaX37Trklc6IY2k8GhveTgUmzK1K4tgr6IzAwe/V9umMgaHVdhkPdk7XUgeH1\nZGDWerJ2gnVGfBk0kjWyDBRDlwEYvgzApFVOy2C1NS7YtVg9WXiVL/xpU8VsbE6E4f1kiMP7ycCk\nnBQpXFsFwxHmFvdGKqihT3dCkC5g1605BEPryTA6rycDs9aTNXPDjPgyAON4QlX9fI4ykDiOMhAm\nrXJaBsutAckarScL/5BbjVy8n0wYh59MGIefDEyY5s/zSTwDUeuB2astnPFIUqADjt2aM2BoPRlG\n5/VkYJZ6smZql8zIUQbCOMqgMp4yqIynDIRZlsF6a7QZOkcnAzfKwNkaN7rZZ6meDJyt4cnW6Xuj\nDJytcaObfZbqycDZGp5snb43ysDZGje62WepngycreHJ1ul7owzErXEUbTe64WepbAaOoo3N1Ol3\nswycJ1Q3u+FnuWwGztZgM3X63SwDq62hXFa0nkwx4Q0qj3dGjQEGrd1tUD0/3t8fb393QLiuGdZP\npuxntJ+smRCnJ1NhaD2ZZlg9mWJiwqZ6siabWA9dBk0csgzAqBvVTMQ4wNzQMro1p9Az3dLZfym9\n2hrKZUXryRQjNrBmXsYBGLSMbs0p1TO9W/yV+E+YwfB+MmU/o/1kakJ/frhPcyEMrycDw+vJwKRs\nch+awnroMlBx6DIAgxvV3HPrAHNDy+qnz0nPdRmstoZyWX0+Xon/5D3Eb5hwHA1iuwcYtHiG95Nh\ndN5PhtFpPxkQXk8GhteTgeH1ZGBigud6Mp1+xdBl0DBhMKYMwOBG6XlYbTBoWf30OfRcl8FqayiX\nlXxiS8ew2mBgA7P66XNg0NLXrTZ6PtLmY/xkYBx+Mhmd9pM1EyI/sgLGoSeTqTn0ZMKEpC70ZDrl\nmNvFl4HEcZSBMLhReh5WG3NDy+qnz6HnugxWWyONlz1bdE7A5M9n6jkt29XnFSwr78uO6mLqyfvJ\nEhkZh58Mo9N+MiAhILk1wPB6MjBpYZSerGEWejKV5EsxdBmA4csAjNwoPQ2zDQYts6M6KT03ZbDd\nGmkrX6yeLM8gMbCBqWnNmzlOvI7WvHe+Uub2fH6J82iH5NGzPYDyk1X7mcNPVpE4F3JrXJoJGKUn\naxlOT6aYlZ6sSSMYvgyEcZSBMCk6VwZg0GombxxIz3UZ7LZGqTlST5bnkRhlAzNmN5xCbaM1dOpO\n1J6knyzRmfH4ycroHgmTmhC7NS7F0HoyMLyeTJilnqzNtDCOMiiMqwwkTghfb247k/EIDFpjr/ZM\n7bkug83WqLKKNPReT5a6ZUbbwNqJWUeIg5bVT5+rPVk/WWQrw/vJ6uiOrVGRGJDdGpoRuU0cYPFo\nGU5PBiY+oWcsdeED8W8vn89veXmcKoPKeMqgMnHJ9UYtlp8ugUGLZ5ZlELfGTNEWQsCzFQ7enrOX\ngJvJFEbbwJrr1gHioGX10+dKT95PFmAZ/ev9/eXx3L/iK6PzfjJB4mTJrdEwpJ6sZTg9GZi1nkzn\nGUw+y5QBGL4MwMS8feopTNtg0Jp2Lhd0z1UZrBRtYajywvdrmugjGKmIR2F6G9iSxAvsaC2BcLHG\niXaf64UoczCRuMLz+4Tmg8lX2M9oPxmQMCa5NVqG+2sXDK8nE2ajJ9PZEIa31F1g+DIAE6KTZQAG\nLT13q933nJXB+glVdWbl9b0RpRR/FeY7KwYxa3rducqA7joYh4Wpv3kpOwjixAEJP5kaPd+svZ9M\nISEGtzUahtSTKYbWk4HZ6MlUvsHwljrF0GWgGLoMwKClpm42h56zMlhvjermSv9+eP2mfslVBjYw\nc4bNycp09qymT39QGd5PhtF5PxlGp/1kQMKUua1xKYbWk4Hh9WRgUjqnejKdbDB8GYDhywAMbpSe\nh9UGg5bVT59Dz3UZrLeGOLNiUig9WVBz5mceraxNT81oCwPa6NWeqozDTyaj834yNTrrJ1MIvTXA\n8HoyMLyeDExK5lRPplMNJr0hhSoDMJDC6TGttmLkRln99DkwaOnrVhs912Ww3hrisuL1ZPB5wQZm\nTbA5J3EU3XQwDoRx+MmE4f1kGJ32kwHh9WTCOPRkwjj0ZGCWerIm32D4MgDDl4Fi5EY1EzEOwKBl\ndGtOoeeyDNZboxnxHJwM3CkDZ2vc6W6ftToycLaGI1mn650ycLbGne72WasjA2drOJJ1ut4pA2dr\n3Olun7U6MnC2hiNZp+udMhC3xlG03emOn7WSGTiKNjJRp9vdMnCeUN3tjp/1khlYbg1oemgH0wUm\nzICVMIkVyiFhAsNLmISJ2SElTMLQEiadAtrBVMN4HEyVgcFpe9cxN7TWULmHqftL7IrWFByY0HNX\nCwOzr4UR2ZfCwMRFmKWw2hpK00M7mBRzsRImMLyECQwvYQKT7in1kQDFsBImhdAOJjC8gwkMDE5p\nYYsvYNBadA/WkSrSQgWgNSENBuPQzLYWxjD7UhiZNB+zFFZbA5oe2sF0geElTGB4CRMYXsIEJuaD\nkzAphpUwNUgIxDiYwPAOJjAwK6UbvfgCBq1Fd9xDVABaNmgxOMczm1rAkJjQrhQsJk7ILoXV1oCm\nh3YwXWBiSOqzT4rhJUyIw0uYwIS5kRImxXCfb3gjJwAAIABJREFUvbvUcngHE8LwDibFJPsGo+IC\ng1Ysjvmj3ENUAFpTaGBCz10tDMy+FgaEKIWBCVOblMJiayhND+1gUgyRjpxbxdASJsWkQRgJU8tw\nEibNkFtDIbSDSTFpOYyKSzFiVsoJnX8Fg9a8d7pSagkVgNaUHJjQs5zjmX0tWGHi54W+pkFkGs0i\nJqWw2Brl847B1uRwMDVqn106yhIaJpxjJEwdQ0mYNMNKmBTDSpiA8A4mMCknlINJMWJWWhRFugQG\nrTWT7yEqAK051zOx564WLCZgq1qwkXUpGMysFFZbQzQ9DtGMMEw6Yp/waBhSwqQZVsKkGFrCpBhW\nwiSIw8EkTE7I6ude6tGnrZqV5OKkgThoTbqW07mWUAFozbmeiT3zOR+zqQUjzLYURmZaCqutcVVN\nD5OOuujKMOmwGFYmoeOE3fWdXlesA06+g+ElTGDSoIyEqSIeB1NlYpD6wf7JKuQ0GN7BBAYtGdBq\n9LXE1ELPxHHzOStCPmcxm1qwkXUpjMy0FJZb4yqaHt7BFJap1D67dEiiFENKmJo4YZzvrFCREe1G\njcNLmPo4jIRJUsA7mJowpIMJjDYr2QvH2ZoC0LhmtfI9RAWgZfXO53omnt3VgsVsasFENqUwMPNS\niFtj/R6qqOmhHUwlXUXts0tH6Z2/VYb8azdBkeElTBLHIWESxiFhkqnl5/Nfy78LcwCECS3zVXbd\nT7djCnqzkr5utUuqwyW0rH7xXLmHqAC0ZsjIYBwfs66FYWpEKfTMohT276GKmh7WwVQXXtQ+rq2R\nGVLCVAJFhpYwgeElTGB4CVNm4tRoB5MKQzuYwNQnYJSKK2Ll9jStMtzwrdxDVABaQ996YmDChV0t\nGMyuFgaEKIWeWZTC6gkVND20gykbeV6KAXSXjpJKxCn/Yk/8LzdgeAkTmBSYkjCBoSVMQHgHExje\nwQQGZqVanLPvYNCa9c3nyz1EBaA1BQcm9NzVwshshVwDQpTCwKRFmKWw2hpK08M6mIKr9OP6+Cyv\nsOzSUXILhpcwgeElTGBSYErCBIaWMAHhHUyKoR1MYGBWKhmdfgOD1rRzulDvISoArRk5MvTWuGT0\nfS0MYYhSGJi0BrMUVlsjGoR+suaWdjAphpcwSRyHhEkYh4QJDC9hAkNLmIDwDibF0A4mMDArzaq1\nngeDVr1mfMc9RAWgZQDhlMXgHM3sagFDYkJ/36RizTgWMy2F5dYwhz8nTwZukYGzNW5xm88i/Rk4\nW8Ofs0PcIgNna9ziNp9F+jNwtoY/Z4e4RQbO1rjFbT6L9GfgbA1/zg5xiwzErbH4by5vkYOzyJMB\nIwOb/+bSIM6pk4FbZOA8obrFbT6L9GfgbA1/zg5xiwwstsYgsyIcZQMT3qHyeH8sP4LXM4ykrGcI\nSdmIiOFseqcLE64PNM+A5pm9nqzOCKP/of1kYOKMTD9ZvJAfQw74MggDCM2WgTBEGcjgNczH+/vj\nbf2f24/MrAymW8OQWW0dZQazlZSNzF5SNjJbSdmI7B1lwkAwxuvJwKhWqbbu2xhnqydTY1ba4SdT\ndJjK8pNTdXQla+PLAHH4MhBmWwbG1NIHB15/dflVhyMzL4PJ1jBlVhtH2YQJE5tLyixmJymzmI2k\nzEJ2jjIwaMEGptKtmuhptVRH1bR6bvRkl8XwfjLQcRq2nyxPED3RuugyaJgwIFMGYDZlgI5obarA\nTNu8DCZbI6Ys+9WUkWv9ccSYzIEhJGUDQ0jKRmYrKRsQwlFWGGthcbXmY2BCL5wzEXW99NzryUaG\n95M1M5r4yWSimHtt0WWAOHwZgNmXQZ2QJPixrYIxbfMy2G4NJbOicwKGkJSVBYJJd2UtKRuYvaSs\nRxhH2ZD7bpJSP6oxMLjZqlfbHJi9nmy8x2lIyk/WzGjiJ5MJDnO76DJAHL4MwKQZLMtgmNq+Coa0\nLcpgtzW0kWvvKMuTVQwhKRuYlJO1pGxktpKyASEcZZmJ8xnoNEnrS8+AtnrncxYTrqz0ZMM9ziMx\nfrLYUyLO/GR5ON1TGLYMVBy6DBSTZrAsA1kEpvZ8fj2WEq+BmZfBbmto7dDeUZYDg2EkZT2TUrJe\nX0kF4gRmJykbwhCOsj6PTcA0z/FLz8QeODf2b6/rnn+f+TkFzZB+Mh1x6ieToJhRbbFlgDh8GYBJ\nE1iXQZ2QSvCuCtStKPS8DDxbI8126SjL4VBBjKSsZ2KU6shIEY0vI7OVlI3I3lGWmRh/pI1ZpVM9\nA3pG1NH7nstXjsZ7nIZn/GQ6ztRPJpO11hMuEmWAOHwZgImtTRmMU9tWgZG2aRnstsZg5Fo6yvJk\nFUNIygYm5GQnKRuYvaRsQAhHWWZwu9TC4knz0TOxE86ZiLquer4tn0zYDOUnUzOa+8lkopgRWuEi\nUQYqDl0GignNTRlgQqW1rwIrbZDWxeh4xK1hK9pKOBi5Xj8j9ngu/kllYAhJ2cCEIOsfmFJrmNte\nUmaFCZGWjrLChG4TOlzpHwMTOuBc3zsf47pqbf7aVT3zS4m8n0xmtPCTyUSHOHwZSJzi8Fq66oY4\ngd6UwYDsq0DdCtAhklUGc0VbQfP0vqIjL22Nt3wjJHVNY2AISdnAEJKynqm/eReSsh4p84atrFlI\nPihMOJjQDAPa6J1PDXGCCj/+61WudhsbGN5PFgbM9MJPJkGHOHwZSBzGVTfE2ZdBjxBVUG8kppYX\napXB7glV+BehCL+FbUE4yspkwRCSsoEhJGUDs5eUDQjhKCsM8oiF5YwaXwcGtNE7nxqZrZ5svMe8\nn6ybkeknk7kOc+PLAHH4MgCzL4NhavsqGNM2L4Pt1lDOrJiUNJLkrW/UyYpni5CUjcxeUjYwe0nZ\ngBCOssrUn7LhV3x4E8LH0so+MqD7dNXjgdnryYx7HF+0/PV8r4Ma34c4qY/pJxN6YAhV3cA4ykBl\n61d+wixT6RtDmH0VjGmbl8Fsa1gyq52jzGJ2kjKT2UjKLGYjKbOQnaMMDFqwgfX3KR+jp9WimZ2e\nTGnQEIf3k4GZ+snKVNETLb4MwPBloJh1GaAjWpsqMNM2LYPZ1rDv4jl7MnCbDJytcZtbfRbqy8DZ\nGr58nd63ycDZGre51WehvgycreHL1+l9mwycrXGbW30W6svA2Rq+fJ3et8lA3BpH0Xab230Wymfg\nKNr4XJ2et8rAeUJ1q9t9Fstn4GwNPlen560ysNgaRWYFW9hHlFktbWtVlQUmJBNSLDuzQ5yGZpmt\nn2wIE0fm9GSYEK8nA8Nb6oQh9GQ1rcIQljqDmfnJkPaaN0mwowyECcOxZSAMVobJtK1harylru1p\nlsF0a1SZlbKF0X4yxYhxq10UjsY4mkY/3RqZrZ9sRNKAyw/LGAytJ1OLoPVkYLZ6MkkrmK2lzmDm\nfrKa7poDlWC6DBQjseu4/fcxjlpZ3zkfj4gOyDKpn1kGk60B7RVsYbyfDAzGsaeK62DQ4pmNn8wK\nEwfn9GRqQrSerGFCJEZPBmajJ1OuMTC8nwzM3E+WU4+8IcF8GYDBOHnc/iuug8Es+97p2EIAm4hK\nW9vTLoPJ1ghDl3fDwxamZG12ZIPBODNkjIOIPLP3kw3LCYOTejI1oc2nUrFcMLyeDMxeTzamjfeT\nIc7cTyapL3lDgvkyAIO8yLh9Y4iDWfZd6/GANAFrr/b7wITLkzLYbg3Ywng/GZgQuEymnaE6KtfB\noKV6tc2BSZeXfjILIfVkakL01gDD68nApOUs9WQ1rWB4P5kwCz+Z5FvfvpRgvgyam6LHkcFVQ19P\ncWSWqlPbHJB0eVkFNW1Nz0kZbLdGnk2whSnxWjtDdaQnWw1jzTnVtzab67CSoVU7qu82s/STGYhD\nTxZipwl59GSFcenJapy4VtYoUpmtpa6pi7SeuZ8sTiA9dN5igr1lUG+KHqeOrb/r65UJ11dlYCMK\n1uOXtsHMyoDbGtEWBruUEdEIXA1jejIWqa9X5rrQIpmdn2wM49GT1Ql59GSZ8enJapyw6LWerCnz\nkiyHn6zEmfvJJOuSt5JgVxmomyLjyMhtQ64rRmWj7ZyPDKSBOWZaBtzWiH/Bu3ISZlX/6pcFWDMN\n5/T1yoC2IZtZ+slGxKMnayZE6sky49OTIU7VY9gJsNLm8ZPVOFM/mYTVeYsJ9pZBvSl6HBlcNfT1\nytRZqm66aSOAdd/aHplpGVBbI9nCfH6ySwxjejJ1gvq7ui4MaN0TbZNZ+8kGxKUnayZE6skK49KT\nIc5GT6Z/ouS0ufxkEmfmJ5Ncq7ylBDvLQG5KM46MjkZz/TsbpWSW6KZbJiIBdU+0B2ZeBnFrLBVt\nYdT3/McnhGiI1LUQuDCRXomrmuuK2fy5izELs/eT9YhTT5Yn5NOTZYa31MVkyMLxC7TLcTns13P5\n/GQSJwxn+ckkaImDBPNlACaMhvnK0E1jiBMZPcumdzoYkCbg2D+e6ZlFGWwVbbCF5ZsVZW3Th6xf\nGcbk3ISS62DQYpm9n6wP49KT1Ql59GSV4S11SPVeT1bvsTD1CdjCUjcwJbuWn0wSX/KGBPNlACaM\nJvmXkdvGEEdW1vZTRwPSBFQdVbNnFmWwf0IltjCHn0yYMCkyJ+U9G+9/6rs3Ymv2qGNKnL2fbEDS\n2JyeTMI49GTCOPRkwuz1ZJJWYRx+MmHmfjLJfMkbEsyXAZgwWs2/DNw1hjiqILqu9XBAmoC1V/t9\nYNJlswy2W0PZwmg/mWLonIBBq12VOioLRM+vrZ9sQNJwlJ4MYdI7EShLHRheT6YY2lIHhveTgZn7\nySTZJW8qwXQZKIYuAzCYpcylawxTA9z1xOHApEtmGcy2hmivlC2M9pMpRsbB7JqWXAeDVtMTBwZz\nbfxkFkLrydSEaD2ZYmg9mWbWejK4xhRD+8kUM/WTlWxL3lSC6TJQDMbBbdQtXJcbqWape0rbQFRA\n6dY0LGZaBrOt0Yx4Dk4G7peBszXud8/PiqkMnK1Bpel0ul8Gzta43z0/K6YycLYGlabT6X4ZOFvj\nfvf8rJjKwNkaVJpOp/tlIG6No2i7330/K95m4Cjatik6He6ZgfOE6p73/ax6m4GzNbYpOh3umYHF\n1hgEWA4/mXi2HH6yynj8ZJUh/GTDcpJqbPH23lAQSiuW363K68liOWWmHSeeHx5DnL2ebJxbax0b\nYsQTQ5x4UmYZD8bHwPBlEAeroz8e74+l3q+PQ5RBjxBVMKZgWgbTrTEKsC7aT6ZUWbSfDAzvJwOz\n9ZONy+H1ZKlc8scUaD2ZYmg9GZitnqwdM81NJSMNNH6pOUCcrjUiRhy+DNTodBmA2ZaBsZz0oYHX\n8B9Yzx4jMy+DydawBFgX7SdTAqzP+NFbxk8GhveTgdn4yazl8HqymOds8eL1ZGAQO54bH+31HGej\nJ1OusTheZpCMMUY8Y8WJ520/WbwyYegyiCOU0ekyALMpA2s5myowUzAvg8nWiIvKH+dT2qv1xxHT\nonqG95MhDu8nU0wy03z8jpOwH8NyeD1ZGLBYvHg9GZjQKrHtiTXXS5y9nkyNWRgkg48Tek78ZDIG\n5l570mWA0fkyALMvg2FqvKUOYeZlsN0aKUdZe+XJSfj44jN8vIj3k6k4sUn5ycDs/WTIY56aT09W\nLF4uPRnMXzp2mnP3BdcLs9eTqa2BOCXp3eg4HOKES5pGT7RGxlMGZXRXGegZLctgmNq+Coa0LcqA\n2xpZe+XzkyXG6SeDXov3k5W5PZ9fS3MT8hj+Go3PST16smLxcunJlPlLx0bRoSXXFROurvRkuMcN\ngwRidLSMOA2NnmiNjKMM6uieMqhMmsKyDIyp7apgTNu8DIitIdorh5+sMC4/mcQJSVlWuXqGAmbn\nJ5M8VsShJ6sWL4+DqTLxDkvsdLvHL/W6ZjZ6MhlTMXVlY4ByZoyj6Ak1MnwZ1NE9ZVCZNJ11GYxT\nu3ZVMKZtXgbE1giTVNor2k8WGa+frMapeozJ7ZIFxuuZ2frJah6B8HqyavHybI3KxHg6djzuH/W6\nZjZ6MhmzZWoC+wj5eIzT0hY1MqkXVQZ1dE8ZVCZG2ZTBOLVtFRhpm5YBtzWu6sy6wr5c/FNAnWzK\nXmR8fjLRazn8ZJnZ+8mGqV0XqycTi5dDTyZMTEUTO+Wm/VKuN8xGT1bHbBlJYDt+PRri9HTtqL4P\nTL7GlAFG58sATIizKYNhavsqsNI2K4O4NZaKNmiveD8ZGN5PBibkJP8jgrpBXbMkBczeTzYgeUhC\nT6YsXrSeTDEhTondLQKH+XrHbP7aHRgkAwN3rZ5pI3ady2HPXHQZqNHpMlBMiL8pg2Fq+yoot6IN\nEyJZZbBVtEF7ldf3RijawPB+MjDxFylVF5cw9Tfvwk+W8wik3HlCT6YsXvleEZY6xYRAJXYJOX7L\n1xumCt7GzuXMwEgypkiZB+KgxTPl5QuiDNTodBkoZl8GfQqIKhhSUBZulcH2CRW0V7yfDAzvJwPD\n+8nA7P1kpTyB8HqylL4E8noyMKFVYpfbMH7T1/ME89srGEtdHC0xWNkYoJwZ4wjNM3wZpDHTtPgy\nALMvg2E5+ypobkWa2rwMtlsD2qv0TgTKTwaG95MphvaTgdn7yUoegfB6snS7ssWL1pMpxrc1Upy9\nnqy5x4nBylJs64uuJVjJ0KIYvgxUDn7CX9cf34t3cIzrif9Ik7ehNa14bljOvgpGZl4Gs61hyaxo\nP5lWZQUJ2OvibWVWnIv1k6k4Gz+ZFYbXk4WXwX6ez/hsjdeTgUFs+x7r6yXOTk8GRVsYss5N5GZ2\nGJMBbUPG3C6+DNTodBkoZl0G1tQ2VWCmYFoGs61hZ+qcPRm4TQbO1rjNrT4L9WXgbA1fvk7v22Tg\nbI3b3OqzUF8Gztbw5ev0vk0Gzta4za0+C/Vl4GwNX75O79tkIG6No2i7ze0+C+UzcBRtfK5Oz1tl\n4DyhutXtPovlM3C2Bp+r0/NWGVhsjSLAgpGL95OBCclUIi0ztUMc3k+GOHs/2RAmTqYKxMyJydQR\nBq0JYTDhXUGknkxGJ/RkNa3COPxkinl7vBerzWxBcvvqKvgyaFYu40wCyfUaZ18GIyLKvkkQuT3N\n1OwymG4NkVlBzMb7ycC0ji9jwmMch59M4mz9ZGOYNJflh2UMRgIaK0mnDIbXk8noWz0Z0irM1lJn\nMHM/WV2frAer4MsADGLXgbvvY5xtGYzItgowDTW1MBOzDCZbQwmwxMjl8JMJo8bpcpEP1XVhHH4y\nYTZ+MitMnACpJ5MwvKVO9eT1ZBJnoyfTrjFhHH4yYeZ+MuP2hFPRtecqg8Ko/OeBu6/qumRrUwYW\nsqmCNm1lanEmdhlMtkbsXz7OJx+4c/jJhFHjxDlYjyGOw08mcfZ+siFMmEuVjlnTiudGRgLOkJFx\n6Mlk9OgOCL/yF2/mH+M4/GQSZ+4nk/WVHGAVfBmACaPVXMrAXWOIsy+DAdlXQZ1GM7VJGTi2hsNP\nJrnnc3IJ4/CTCZPynEVyXcrrYb01GtEysNpPfx8ZTeueaA+MQ0/Wjr7Uk9V7jLQ5/GQ1zsJPJgsq\n68Eq+DIAE0areZGBu8YQZ18GA5KGXFZBnUYztUkZEFvj6+v9Lfz88vjJtMWLzMmlmbBGyk/WMks/\nWZ2GQhoZWHen0uHIKNoCwrmBcejJ2tGXerIxTvjo/8ZSNzJzP5msrqxHVuEoA2HiYDUvMnDX6OOU\ny6sysJFlFdRp6KnNymC/NaqRyyNhqkxcH5mTSzOsn0wxOz9ZnQaQRgbW3ah8ODKgTUAtt/b06Mkq\nk8Ze68kkrYqh/WTCzP1ksrycA6yCLwMwcbCaSxm4a/Rx8uW/z/zcsuucDw1kVwVlGnpq0zLYb400\njWDk4nNS1lEsXlxOWmbymkHp1OVZbGFLP1kzjYRoGZiM3DRGJl2WgE3nctAzHj2ZHr3qMawY8Vwf\nJ/x9/Pby+fz+OwNMZuonk1FyHKyCLwMwcbBmvjI6Gn2cfMV86ahCNrKsAtka6Q/pt2jnmZYBuTWC\nkcvhJ8tzLxYvLidludX89cY9m0hUZdZ+smYaEWlkYDXb7feByZcRsO2ejgaG15M1o2/0ZG2pxRnx\nfjIdZ+Ynk5WV9cgqHGUgTBysyYuMjkYfJ11Zl4GJrKugTgNTm5dB3BpLRRuMXLSfrLV4kTlBnJjH\n+ociUte0yphg9n6yHhksXU2AfNAz7cIMIJwaGFpP1o6+/Hlpxgk/AaOgZqHtGOYWCdtPlq+Er4XB\nKvgyAINxZNy+McSJzLoMBmRfBcNyFmWwVbTl9UUjV75ZhJ+stXiVBfSpkONmgdn8RfrJVJy9n6wP\n08jAZDJto2dUwLajOhqY/GLI13PxTGdg4ose67Jo73FIW30CtrDUDUyZtuUnkxWVuWEVfBmACaOV\ncWTcvjHECbKd6NrJ5d73TscDsq+COg2Z2qIMtk+oYOTi/WRgwhLInCiG9pOB2fvJyjSApOxmzkw8\npg4GrRlSl4uevJ4MzF5PNsbh/WSIM/eTyfpK3rAKvgzAIJcybt8Y4pR3byw0dQOyr4KatmZqxXHX\nz2i7NZSRi/aTKYbeGmB4PxmYvZ+s5BFIygSlJwODVp9HOR7j0HoyPfpOT1bvMRjeTwZm7ifr16Nc\ne3QZKIYuAzD7MiipBrKvgpo2MPMymG0NCLBg5OL9ZGAwjiS7aeC6MA4/mTBK1tYMXw+MMOFSlZvV\nXu13i0HAtm89spiL1pOp0dd6MuUaA8P7ycBM/WRlQVgPVsGXARg1Tk1V811dr9nalYGB7KpApQ1T\nm5bBbGs0Ez8HJwP3y8DZGve752fFVAbO1qDSdDrdLwNna9zvnp8VUxk4W4NK0+l0vwycrXG/e35W\nTGXgbA0qTafT/TIQt8ZRtN3vvp8VbzNwFG3bFJ0O98zAeUJ1z/t+Vr3NwNka2xSdDvfMwGJrDAIs\nmL2muRqY8A4V0k8mPR1+MmEIP9kwtaQa+zNdSryQGTUhWk+mGBlnHmqIs9eTjXPbW+pGJk6J0tSp\n9dBloBjcqGkS+hw0tE31CFEFYwqmZTDdGqMA64INzJ6pJcBqVVgGN8Zx+Mkwenqj9ev8YzxjGF5P\npiZE68kUg7QYy4+n6tzAbPVkBpPeT7vUaYxx0oSWH5oyGLoMsJ4LN4rOgaJtxpjargqMtM3LYLI1\nLAGWso6Zc50woW80e9kPi3H5ycroGz+ZFYbXk2FCvJ4MjIpt5gDXwWz0ZMo1Bob3k4GJE7L9ZHmq\n1tz4MlBxRLxmpsBcj6ItyJrapgrMMPMymGyNmLKsaFMyq83nzixG0db64rkhDu8nw+h7P9kQhteT\nYUK8ngwMljjLQE0BmL2ebGR4PxnihBlN/GQy15I3xdBlAAY3SsbtG0Mc0H3Xejwi6Ufwx+/awfg+\nMPMy2G6N/FHBFITOCRi0jGnmU2WyXU/KTwZm7yfrwzj0ZGmeaUK8ngxMaJXYeb3GV309xdnryZox\nkazlEyqTmfjJZJbD3CCFkz59Y2Bwo/qucjww6QpWJh3RGJB9FQwpWJTBdmsomVVrEMMU0SqTBYMW\nOnWtgUnXp0/B0tWR2frJBiR/6P3x/Ormow4LIxNy6MmEiQ09TrrQfdHX1cJXerJmTDCUnyxFL8zM\nTyYzHOdGlwHi8GUApmvJhNAwpkZb6uIoKQXzMthtDS2zErMXZte18mTBoNV1VIc9ky5RfrJm9J2f\nbAhD68kwId7BBCa29D1MV7ov6rpa+FJPpsesDOknS8ELM/WTyQzHubFlgDjNjZKR28YYJ1yvK2u7\n1iMD2VXBmLZ5Gey3BmRWaUqEnwxuLrTqcsbveYFtz6rHGHvnMyOz9ZONCKsniyHzhHxbA4tQ99Bc\nEa6D4S11mmH8ZFhPqL1w8IlfOsbs7LldRBkgTntzjSDhlBVHr8ygRmRbBUaYaRnstsYFmVWeHOMn\nA4OWsbR8qiyw6Un6ycDs/WRGGFZPFuaZJ+TQkwkTF4l7GI/GB66rha/1ZGpMxXB+sjCBwsz9ZDJJ\nc27XxZQB4uBGybh9w4rTrKwHdFpLx30VWGmblUHcGktFGzxbEKKNsyxnyvrAoMUzoefypfZwfYyz\n9ZMNSJ7Qa34dzp4dbleZEK8naxaBccg4odtGT1ZTgDi8n0yYhZ9MJoq5lxzwZSBxir/ri/u7Djcf\nLZmObgxTi55OzlIXujWDW2WwVbTlFxeiYSyXeZao6RmqdpksGLRUr7Y5MLyfTEavv3kXfjIjTJwH\noyfDhHI6GUsdmBAE97Bdej2S6xCz7fRkGLMyvJ9M5rbwk9WpjXH4MpA4l9woGXZoGDmoKxv6lhM9\nQlTBuJw8llUG2ydUkFnB7DWbaw0MBi2e4f1kGH3vJyt5BMLryTAhXk8GJixc7uEkCXL9Jf9aCnBK\n9kJPhjErw/vJmrlN/GQy02FufBkgDrIu4/aNIQ7ovms9HpB9FYxpm5fBdmtAZgWzV53b8L1OFlYy\ntIbO5cTIXKyfDHPb+8mGMA49GSbE68nAOLZGXfheT4Z7XBmHn0zP7booTV3QbJZt6ygDYXCjZlUw\nrgcRZ0y9pRJmXwVjmHkZzLaGJcCC2cuerMVoFZZFmQztJ8PoGz+ZFcahJ5MJOfRkwqjYVga0N6wy\nOz2ZxXj8ZDVOmA+rqbuEcZSBMLhRZgrM9SCiyai0SphNFZhhpmUw2xrmbM7Jk4H7ZOBsjfvc67NS\nVwbO1nCl63S+TwbO1rjPvT4rdWXgbA1Xuk7n+2TgbI373OuzUlcGztZwpet0vk8G4taw30N1nxyc\nlZ4MGBmYv4fK6HxOnQzcJwPnCdV97vVZqSsDZ2u40nU63ycDi62RBVhKe0X7yRQTMlnGmeZ0iEP7\nyVScvZ9sCBMnROnJVBhaT6YYWk8GhtCT1bR+vL8/3tL/SY4Wn+qpnwxD1Ly9Pd6z/8ZRBsKE4dgy\nEGZfBsPU9lVQp9Eu3CyD6daoAqwL2ivaT6YYkWIh121rjMP7yRBn6ycbw6RpNJ9naScGdRrCOCx1\nSButJ0OcrZ4MaUUctPqFlOMxB3M/WR1A9vWrAAAP6UlEQVSiMqonXQaKwXzrwN33Mc62DEZkWwUy\nDT21MBOzDCZbAwIsaK94PxkYjNOlohziOhjeTwZm4yezwsQZcHoyhOH1ZA0TIjGWOjAbPZlyjYFB\nq+S2+2blYO4nyzAY9OTLAAzG6SZVDnEdzKYMLGRTBSptCBNnYJfBZGvE/vmjoZCfOfxk6YP4RZUl\n76ovaei/DXEcfjKJs/eTDWHCPFg9mYThHUxIm0NPJnH2erLx9iBin2E5HnIw95P1DHryZQAmjEaW\nAZh9GZQhgeyroE4DTJjapAy2WwPaK95PBobPCRjeTwYm3cmln6zksUFIPZliaEsdmPzhzzTB2Rdr\nbte11JPVe4w4aM3CDMzCTyZj5LmpnnQZKCaMVtYo4/aNIc6+DAYkjbmsgjKNdmqTMthujavKzzx+\nssrEuXI5uSROyRnlJ9Nxwh966dl24ftvdRoKofVkYHg9mTAOPZkwafLpd0i/DDke19PS0lE1BiZ/\nzpvR1ElPRxkIE6dQY6vpNM1yvWFCh1UZ2MiyCuo0dJhZGey3RvCqPONLAS4Jk1JlkTmJ/pYUJyeM\n9JOB2fnJZBqCOPRkwjj0ZIVx6ckkTkjCWk9W77FOm6absqsHQw7mfrKK1DjS01EGwsTBJLaM3DbK\n9YYJvzmf+bll27ccGciuCuo0VJhpGey3RtVeOXJyVSaugcxJw0xeM5AE1TF1nFBM3y/SY2iMCK8n\na8NwerLKePRklYmTr3qMYSHlxLgeTdvUyEz9ZDJAZWpPTxlUJg5Wx5GBu0a9rplNGdjIsgpkGggz\nLYO4Nez/5rIEFu2Vw08mjCMnDXNxfrKWWfvJ+uVcvJ6sD8PoycDwejIwsZQ+u9rpDof1NHTXuRwO\nzHXN/GQyQGGkp6MMhImDyTgyctuQ63pG6zIwkXUVYBo1zLwM5v/NZQkM7RXvJwPD56RlNn/uDnPb\n+8l6xKEnw9R4PRmY/LSW0ZOBCWkzX2pX1dSv52po1VE1ByZfs/xkQhUGPfkyABNazTgyOhrN9TKj\njaZuQPZV0E4jhFmUwe4JVf29HuRntJ9MMXROGob0kylm7yfLeQTC68nA8JY6xdB6MsVAboba6Vr9\nehq661sPe6act/xkFWlrKfaky6AdvaljGR2N5nqe0a4MBmRfBcNyFmWw2xoXtFe8nwxMWHqzAKRC\nWuW6Ymg/GZi9n2wME2eQOZlL1xiY/BLY79VznYHh9WRYzl5PVtMKBq1uFTgc5jb3kwlUGPTkywBM\nGK2MI+P2jSHOXlM3IPsqqNNopjYpg+3WUNor2k+mGDonYHg/GZi9n6zkEUi6N5SeDAyvJwPD68kU\nI9axvoLkeFiPpqVX2xiYuZ9MwMKonnQZKIYuAzD7Mhimtq+CujUQJi3ULIPZ1oAAC9or3k8GBuNI\nspsGrgvj8JMJs/OTGWHCNFg9GcLwejIwvJ5MMWIda7IlB9Z6FC39dMNipn6yAoJBT74MwGAcPSO0\ncV2YXRkYyK4KlKJNwoQ5TMpgtjUw69M6GbhlBs7WuOVtP4veZ+BsjX2OTo9bZuBsjVve9rPofQbO\n1tjn6PS4ZQbO1rjlbT+L3mfgbI19jk6PW2Ygbo2jaLvlrT+LXmfgKNrW+TlXb5uB84Tqtrf+LHyd\ngbM11vk5V2+bgcXWGARYDj+ZeLbCO1Qe74/Fp+/E3SXSLI+fTOLs/WTDciTg5OZDEQYpGVo2ZDGh\nZ4nNMxjHwTwej59fdv94FmO2qzD9ZGUYi9mVgcXsysBgdmVgIHDcTbJgMlJDHTTdGqMAi/eTKQEW\n7ScDw/vJwMBu1q2vHo7LUXDt1HxXijBIydBq+tYDkxEvWO3VfrcYda7tXI7UdZnRzk9mMWm4xYem\nTOb7GR7v5rziSZPZlIHFbMrAQnZVYDHzMphsDUuAxfvJlADrM370lvGTgeH9ZGA2fjJrOYDtuwxF\nGKRkaPEMYvMMYvPMzk+GMdtV2H6yHNdkPh+v73/secWzEyZcmZeBxWzKwEI2VWBObV4Gk60RFlLe\nDa9kVptPpVoM7ydDnOyQWD4JG+a295MNCAKGiRsPKMIgJUPLAMIpiwmnS2yewTg8s/OTYcxmFRM/\nWY5rMpsysJhdGVjMpgxMRBx3fNrmZbDbGlpmtclJmE6uAcXQfjLFpGUxfjLF7P1k/dQUbOcRijBI\nydDiGUmLjVxWHJyzIVzvZrTwk02YiZ8sxzWZTRlYzK4MLCbNYF4GFrKrAoNZlMFua5SPQyeLF+0n\ngwCL95OBSTmh/GSK2frJyk9uIGjZ5ZfPBkUYpGRorZCoFWt7Ln9rGHHK6Cs9WewyxNlY6ixm5icr\nU0jfujj7MhjmRpTBwKTQ6zIYUrCtgjhou5x5GWy3hpJZ0X4yYRx+MmFSSjg/mWZ2frJSnkDQShHt\nL1ERBvMSWnbvfLZl4rn91hiZjZ4sjDowWz/ZyEz9ZHkp6WsfZ18GfRymDHomhd6UwZiCXRWMaZuX\nwXZrXJBZpdn+/E7f7C+1Birj8ZNVJo5c9Rh2FNQamK2frJ9a+AkbXmX5tZRaZn0GNgRas4nF8/EV\nn7ZnjT2nRiaPMyesOOEl0pWlzpjb1E+mAltzW2vq+rkxZdAzcQK7MhhSva0CIwXTMohbY6loazxb\nYWTGTwaG95OBCUFIPxmYvZ9MyrO6uQCHgJNHUoRBSobWpH883THxlMSOB9bDYHaWOitOGPo766is\nIOFcH2fuJ8MAPZOvLMtgiEOUwcCEOLsy6Ke2r4IxBfMySIq239Zrcc39jM4s3k+WsxeZ/EyO8ZOB\nCa3FS+2p3zC3vZ9sQJqA+aD/WhRhkJKh1XeV44EJV5rY0hMNk9n8tTswjJ8sjymrWPjJZHJDHKYM\nujhUGfRMmMGmDIap7avgGpi80OKEk1XHxq/w/Ojlaf1rdXM/ozMrl/lb/m83mkHkYGBoP1kZIru5\nPvLrLjLs0Ojj1N+8QSQ3e/RI6ZcD2lBVhOU79BUWjpZNhHxGVVWo0qZnE3skLaaeG3vnM/U64uz9\nZAOz8JPVuANDlMHI7MtgZOLzqbxd6ly67wNCVMHAlDGtMtj+rQGZFe8nA8P7ycDwfjIwez9ZKU8g\naHUpr4diioOUDK3aqftuMKHHemtYjJzrxq+Hch0z2vrJDCYNt9LUGcy2DAxmWwYGsysDA9lWgcHM\ny2C7NSCz4v1kYHg/mWJoPxmYvZ+slCcQtGq9td+VIgxSMrTazuXIZDZbw2LUOTrOzk+mxmxXYfrJ\ncliL2ZWBxezKwGTWZWAhuyqwmHkZzLaGJcDi/WRKgBWar9bztXLLrTgX7SdDnI2fzAoD2Cw/pQiD\nlAwtnkFsmlGxaWbnJ1Nj6lVM/GQ5rMlsysBkNpo6m1mWgYlsqsBkpmUw2xrmHTknTwbuk4GzNe5z\nr89KXRk4W8OVrtP5Phk4W+M+9/qs1JWBszVc6Tqd75OBszXuc6/PSl0ZOFvDla7T+T4ZOFvjPvf6\nrNSVgbM1XOk6ne+TgbM17nOvz0pdGThbw5Wu0/k+GZhtDcisQi603WzxfiiLAW3n1GKacwbWXC9z\n+7P2k1lIHDm/FdOIEU5ZDK8nCwOImG1tqTPi7PRk5tx2ljojzsVr6kC3grcxd+iJls7GSOhcg0HL\nIlQK0HFTBSYTBzfLYLI1lMwKhjGxgZlT1W4uMKplUVYcfY5lNn4yPWSVtaWh88cdrCj2cng9GRbO\n68mE2ejJ7Lmld4u/TvWFVg7mfrKcEjBoXXQZKEZWZqZarQcMWiaDy2htqsAMkwY3y2CyNSDAgmGs\nNXuN87UY0GP/eMZicI5nNn4yDNlOiNOTKYbWkzVMWAejJwOz0ZOZaeP9ZIgz95Pl1CNvaPFlAAYR\nd7cUDFomg8tobapApQ1MHNwug8nWgAArkvmDfY3Zy5iuxYA2gHDKYppzBtZcL3Pb+MksJIxM6slC\nz/phpOXHzuzl8HoyxNnoycy07Sx1Vg7mfrKcdzBo8WUABivL4/Zf0dNq9b3TsdVxUwUqbaDDYJMy\nmGwNyKwCWuqis4ENE7YY0EP3dMJimnMG1lyvNRv7zf1kE4TUk4Wha5jN1rDi8HoyFScuZ64nU1o3\nMLyfTJiFnyxOIPy1lH4ovoYPSKPFlwGYMFTNYBq3/4KeVqvvnY6nHedVoBYBOgw2KYPJ1sizKY6w\nvKrWOmZON55smHhimZPYwWBwLl23vgxxwo3Mn820eqdzPULqySJbF8HqyRRD68kUE5vXTk/Wp5r1\nk6k4cz9ZmkH5UvIWjkLLVwaZiePUDMb25NHEKX1wzoRwWVq7KqhpC+NlZlYGq60R3VzxkVfVqpXS\nBeNLy4A2usqpkcm+LulgNAZm7yfrEVZPFqPXG8vqycDwejIwsXXt9GTJT9YwpJ9MMXM/WexUHjVv\npkStduq+t0y8WDPYdVSHI7MtgwHZV4GkrQ4+LYPV1qh/t3u2RstwORmZsKE3PzPrdZ3xjZ+sR1g9\n2bCIpaUOU89T4/VkbZyqx1C10zX79bB+Mh1n6idTsWqcvDLuJ2TL6Ihq4K45Mshl17UeWsimClRl\nZXpaBoutkQRYcRL1HqffIX8X/7VCEWCBaVvxyHh0cVIPOWf0j6fkut4aaz9Zj9B6shivCUPoyRRD\n68kUEyN+xuPFo18P7SeLY8p6IKybhZI4lkRtAnVME3GC4JaaLZMywoR+hKUujZbpeRnMt0aRWYVR\nSh7F7GVOM50cGNBTyGQ2f+4OzN5P1iMOPRkWwevJwDCWOlmulKz64WZnrl/PxfvJMLc8tOUnk6Aq\nTp6lrwyMlcnQTWOIE67iXNO1HuByae2rQA2ZmUUZTLdGlVnFGeYXb/Pvn2gqmz1GBrSHUeOYmLpe\nt22a1c9zqmgbEIeeDIvYW+qGOBevJ0OcrZ5MVHDC1CdgC0vdOLecXctPlq9AOYeWpwxmEevo+I6e\nVgv9VGvsyFvqZDmLMphtDZFZhbmU8oMNTM1PNw0GtO6o2hajz6mu0tTXy9x2fjIDScNRerLQs4Th\n9WRgeD0ZmJ2erLyzIf/nSDUF2c76Y0la01qNHMz9ZImo76CIcYR2lIEwemVl5PYbelqttm85Mjru\nqgCLUK04nFkGk62hZFZSF1dr9hrmazKgh/7xhMU05wyquV7qYuMns5A0Mqcnk0U49GTCOPRkYNZ6\nMjNtDj+ZxJn7yXLekTe0+DJQjEQ07mc4hZ5Wy2SsjpsqMMOkwc0ymGwNJbOCYUybvYzpmgxog4j/\nppQe8T8mkJ7qHM1s/GRqSAkTxmb1ZGB4PRkYXk+mmKWezEzbxfvJEGfqJ8uZR97QuugyUAwibm4p\nGLRoZFMFKm3N4JMymGwNczbn5MnAjTJwtsaNbvZZqicDZ2t4snX63igDZ2vc6GafpXoycLaGJ1un\n740ycLbGjW72WaonA3lrpFdQ5d/0PfzpezLw72Xgd/43hfg/wKbH339viWdFJwP/JQPhTZfxcf0f\n6Boj4ZsTm9kAAAAASUVORK5CYII=\n",
      "text/latex": [
       "$$\\left[\\begin{array}{ccccccccccccccc}0 & 1 & 2 & 3 & 4 & 5 & 6 & 7 & 8 & 9 & 10 & 11 & 12 & 13 & 14\\\\15 & 16 & 17 & 18 & 19 & 20 & 21 & 22 & 23 & 24 & 25 & 26 & 27 & 28 & 29\\\\30 & 31 & 32 & 33 & 34 & 35 & 36 & 37 & 38 & 39 & 40 & 41 & 42 & 43 & 44\\\\45 & 46 & 47 & 48 & 49 & 50 & 51 & 52 & 53 & 54 & 55 & 56 & 57 & 58 & 59\\\\60 & 61 & 62 & 63 & 64 & 65 & 66 & 67 & 68 & 69 & 70 & 71 & 72 & 73 & 74\\\\75 & 76 & 77 & 78 & 79 & 80 & 81 & 82 & 83 & 84 & 85 & 86 & 87 & 88 & 89\\\\90 & 91 & 92 & 93 & 94 & 95 & 96 & 97 & 98 & 99 & 100 & 101 & 102 & 103 & 104\\\\105 & 106 & 107 & 108 & 109 & 110 & 111 & 112 & 113 & 114 & 115 & 116 & 117 & 118 & 119\\\\120 & 121 & 122 & 123 & 124 & 125 & 126 & 127 & 128 & 129 & 130 & 131 & 132 & 133 & 134\\\\135 & 136 & 137 & 138 & 139 & 140 & 141 & 142 & 143 & 144 & 145 & 146 & 147 & 148 & 149\\\\150 & 151 & 152 & 153 & 154 & 155 & 156 & 157 & 158 & 159 & 160 & 161 & 162 & 163 & 164\\\\165 & 166 & 167 & 168 & 169 & 170 & 171 & 172 & 173 & 174 & 175 & 176 & 177 & 178 & 179\\\\180 & 181 & 182 & 183 & 184 & 185 & 186 & 187 & 188 & 189 & 190 & 191 & 192 & 193 & 194\\\\195 & 196 & 197 & 198 & 199 & 200 & 201 & 202 & 203 & 204 & 205 & 206 & 207 & 208 & 209\\\\210 & 211 & 212 & 213 & 214 & 215 & 216 & 217 & 218 & 219 & 220 & 221 & 222 & 223 & 224\\end{array}\\right]$$"
      ],
      "text/plain": [
       "⎡ 0    1    2    3    4    5    6    7    8    9   10   11   12   13   14 ⎤\n",
       "⎢                                                                         ⎥\n",
       "⎢15   16   17   18   19   20   21   22   23   24   25   26   27   28   29 ⎥\n",
       "⎢                                                                         ⎥\n",
       "⎢30   31   32   33   34   35   36   37   38   39   40   41   42   43   44 ⎥\n",
       "⎢                                                                         ⎥\n",
       "⎢45   46   47   48   49   50   51   52   53   54   55   56   57   58   59 ⎥\n",
       "⎢                                                                         ⎥\n",
       "⎢60   61   62   63   64   65   66   67   68   69   70   71   72   73   74 ⎥\n",
       "⎢                                                                         ⎥\n",
       "⎢75   76   77   78   79   80   81   82   83   84   85   86   87   88   89 ⎥\n",
       "⎢                                                                         ⎥\n",
       "⎢90   91   92   93   94   95   96   97   98   99   100  101  102  103  104⎥\n",
       "⎢                                                                         ⎥\n",
       "⎢105  106  107  108  109  110  111  112  113  114  115  116  117  118  119⎥\n",
       "⎢                                                                         ⎥\n",
       "⎢120  121  122  123  124  125  126  127  128  129  130  131  132  133  134⎥\n",
       "⎢                                                                         ⎥\n",
       "⎢135  136  137  138  139  140  141  142  143  144  145  146  147  148  149⎥\n",
       "⎢                                                                         ⎥\n",
       "⎢150  151  152  153  154  155  156  157  158  159  160  161  162  163  164⎥\n",
       "⎢                                                                         ⎥\n",
       "⎢165  166  167  168  169  170  171  172  173  174  175  176  177  178  179⎥\n",
       "⎢                                                                         ⎥\n",
       "⎢180  181  182  183  184  185  186  187  188  189  190  191  192  193  194⎥\n",
       "⎢                                                                         ⎥\n",
       "⎢195  196  197  198  199  200  201  202  203  204  205  206  207  208  209⎥\n",
       "⎢                                                                         ⎥\n",
       "⎣210  211  212  213  214  215  216  217  218  219  220  221  222  223  224⎦"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "globInd = []\n",
    "for rows in range(len(xVec)):\n",
    "    cols = []\n",
    "    for col in range(rows*len(xVec), (rows+1)*len(xVec)):\n",
    "        cols.append(col)\n",
    "    globInd.append(cols)\n",
    "\n",
    "display(Matrix(globInd))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Make the inversion matrix"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Make the diagonal for the inner points"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "c = []\n",
    "for x in range(nx):\n",
    "    indexStart = (startXIndex+1)+(startXIndex*x)     # Multiply by 10 due to index system\n",
    "    indexEnd   = (startXIndex+nz+1)+(startXIndex*x)  # Multiply by 10 due to index system\n",
    "    for ind in range(indexStart, indexEnd):\n",
    "        c.append(symbols('c_'+str(ind)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Make the ghost points"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# The inner ghost\n",
    "innerGhostStart = startZIndex\n",
    "innerGhostEnd   = nz\n",
    "ig = []\n",
    "# +1 in the range as last point is not included\n",
    "for z in range(innerGhostStart, innerGhostEnd+1):\n",
    "    ig.append(symbols('ig_0_'+str(z)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# The outer ghost\n",
    "# nx+1 as we want to go past the last inner x grid point\n",
    "outerGhostStart = startXIndex*(nx+1) + startZIndex\n",
    "outerGhostEnd   = startXIndex*(nx+1) + nz\n",
    "og = []\n",
    "# +1 in the range as last point is not included\n",
    "for z in range(outerGhostStart, outerGhostEnd+1):\n",
    "    og.append(symbols('og_'+str(z)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Create the diagonal of the matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "InvM = diag(*ig, *c, *og)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Create the non-diagonal values"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Set z+1 values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "for x in range(nx):\n",
    "    # The indices referring to the matrix index\n",
    "    # The last -1 is there as the matrix indices count from 0\n",
    "    startRow = (nz+1)+(x*nz)-1  # Starting at row+1 after inner ghost point sub-matrix\n",
    "    endRow   = (nz+1)+(x*nz)+(nz-1)-1 # Ending row-1 before the last z-index (last z will be wrapped around)\n",
    "    # +1 in range as last point is not included\n",
    "    rows     = range(startRow,   endRow+1)\n",
    "    cols     = range(startRow+1, endRow+1) # Column is shifted +1 from the diagonal\n",
    "    \n",
    "    # The indices referring to the spatial point in the grid\n",
    "    # The last \"+1\" is fue to the fact that the column is shifted +1 from the diagonal\n",
    "    startInd = (startXIndex+startZIndex) + (startXIndex*x) + 1\n",
    "    endInd   = (startXIndex+startZIndex) + (nz-1) + (startXIndex*x) + 1 # Wrap around last point\n",
    "    # +1 in range as last point is not included\n",
    "    inds     = range(startInd, endInd+1)\n",
    "    \n",
    "    for rInd, cInd, ind in zip(rows, cols, inds):\n",
    "        InvM[rInd, cInd] = symbols('zp_'+str(ind))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# The wrap around\n",
    "# The index referring to the spatial point in the grid\n",
    "startInd = startXIndex+startZIndex\n",
    "# The indices referring to the matrix index\n",
    "# Last -1 as the matrix indices are counted from 0 \n",
    "startRow = (nz+1) + (nz-1) - 1 # nz+1 below from the ghost sub matrix, nz-1 below after that\n",
    "startCol = (nz+1)-1            # nz+1 left of the ghost sub matrix\n",
    "for wrap in range(nx):\n",
    "    row = startRow+wrap*nz\n",
    "    col = startCol+wrap*nz\n",
    "    InvM[row, col] = symbols('zp_'+str(startInd+startXIndex*wrap))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Set z-1 values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "for x in range(nx):\n",
    "    # The indices referring to the matrix index\n",
    "    # The last -1 is there as the matrix indices count from 0\n",
    "    startRow = (nz+1)+(x*nz)-1  # Starting at row+1 after inner ghost point sub-matrix\n",
    "    endRow   = (nz+1)+(x*nz)+(nz-1)-1 # Ending row-1 before the last z-index (last z will be wrapped around)\n",
    "    # +1 in range as last point is not included\n",
    "    rows     = range(startRow+1, endRow+1) # Row is shifted +1 from the diagonal\n",
    "    cols     = range(startRow,   endRow+1)\n",
    "    \n",
    "    # The indices referring to the spatial point in the grid\n",
    "    startInd = (startXIndex+startZIndex) + (startXIndex*x)\n",
    "    endInd   = (startXIndex+startZIndex) + (nz-1) + (startXIndex*x) # Wrap around last point\n",
    "    # +1 in range as last point is not included\n",
    "    inds     = range(startInd, endInd+1)\n",
    "    \n",
    "    for rInd, cInd, ind in zip(rows, cols, inds):\n",
    "        InvM[rInd, cInd] = symbols('zm_'+str(ind))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# The wrap around\n",
    "# The index referring to the spatial point in the grid\n",
    "startInd = startXIndex+startZIndex+(nz-1) # +(nz-1) as this will be the last z point for the current x\n",
    "# The indices referring to the matrix index\n",
    "# Last -1 as the matrix indices are counted from 0 \n",
    "startRow = (nz+1)-1            # nz+1 below the ghost sub matrix\n",
    "startCol = (nz+1) + (nz-1) - 1 # nz+1 left from the ghost sub matrix, nz-1 left after that\n",
    "for wrap in range(nx):\n",
    "    row = startRow+wrap*nz\n",
    "    col = startCol+wrap*nz\n",
    "    InvM[row, col] = symbols('zm_'+str(startInd+startXIndex*wrap))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Set x+1 values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Indices referring to the spatial points in the grid\n",
    "startInd = startXIndex*2 + startZIndex # *2 as we start at the second inner x-index\n",
    "endInd   = startInd + (startZIndex*nz) # *nz as this is the last z-index in the current x-index\n",
    "\n",
    "for x in range(nx):\n",
    "    # The indices referring to the matrix index\n",
    "    # The last -1 as the matrix indices counts from 0\n",
    "    startRow = (nz+1)+(x*nz)-1      # Starting at row+1 after inner ghost point sub-matrix\n",
    "    endRow   = (nz+1)+(x*nz)+(nz)-1 # Ending at the row referring to the last z-index\n",
    "    # Not +1 in range as we do not want to include last point\n",
    "    rows     = range(startRow,    endRow)\n",
    "    cols     = range(startRow+nz, endRow+nz)           # Start at first index after last z-index\n",
    "    \n",
    "    # Indices referring to the spatial points in the grid\n",
    "    inds     = range(startInd+startXIndex*x, endInd+startXIndex*x)\n",
    "    \n",
    "    for rInd, cInd, ind in zip(rows, cols, inds):\n",
    "        InvM[rInd, cInd] = symbols('xp_'+str(ind))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# x+1 for inner ghost point\n",
    "# Indices referring to the spatial points in the grid\n",
    "startInd = startXIndex + startZIndex   # First inner point for first z\n",
    "endInd   = startInd + (startZIndex*nz) # First inner point for last z\n",
    "\n",
    "# The indices referring to the matrix index\n",
    "# The last -1 as the matrix indices counts from 0\n",
    "startRow = startZIndex-1    # Starting at first row\n",
    "endRow   = startZIndex+nz-1 # Ending at the row referring to the last z-index\n",
    "# Not +1 in range as we do not want to include last point\n",
    "rows     = range(startRow,    endRow)\n",
    "cols     = range(startRow+nz, endRow+nz)           # Start at first index after last z-index\n",
    "    \n",
    "# Indices referring to the spatial points in the grid\n",
    "inds     = range(startInd, endInd)\n",
    "    \n",
    "for rInd, cInd, ind in zip(rows, cols, inds):\n",
    "    InvM[rInd, cInd] = symbols('igxp_'+str(ind))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Set x-1 values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Indices referring to the spatial points in the grid\n",
    "startInd = startZIndex\n",
    "endInd   = startInd + nz\n",
    "\n",
    "for x in range(nx):\n",
    "    # The indices referring to the matrix index\n",
    "    # Note that x starts counting from zero, so we must add 1 to x in the rows\n",
    "    startRow = ((x+1)*nz)      # Starting at row+1 after inner ghost point sub-matrix\n",
    "    endRow   = ((x+1)*nz)+(nz) # Ending at the row referring to the last z-index\n",
    "    # Not +1 in range as we do not want to include last point\n",
    "    rows = range(startRow,    endRow)\n",
    "    cols = range(startRow-nz, endRow-nz)           # Start at first index after last z-index\n",
    "    \n",
    "    # Indices referring to the spatial points in the grid\n",
    "    inds = range(startInd+startXIndex*x, endInd+startXIndex*x)\n",
    "    \n",
    "    for rInd, cInd, ind in zip(rows, cols, inds):\n",
    "        if (ind) < startXIndex:\n",
    "            ind = '0'+str(ind)\n",
    "\n",
    "        InvM[rInd, cInd] = symbols('xm_'+str(ind))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# x-1 for inner ghost point\n",
    "# Indices referring to the spatial points in the grid\n",
    "startInd = startXIndex*nx + startZIndex   # Last inner point for first z\n",
    "endInd   = startInd + (startZIndex*nz)    # Last inner point for last z\n",
    "\n",
    "# The indices referring to the matrix index\n",
    "# The last -1 as the matrix indices counts from 0\n",
    "startRow = len(xVec)-nz-1    # Starting at last inner point row\n",
    "endRow   = len(xVec)-1       # Ending at the last row\n",
    "# +1 in range as last point is not included\n",
    "rows     = range(startRow+1,    endRow+1)\n",
    "cols     = range(startRow-nz+1, endRow-nz+1) # Start at first index after last z-index\n",
    "    \n",
    "# Indices referring to the spatial points in the grid\n",
    "inds     = range(startInd, endInd)\n",
    "    \n",
    "for rInd, cInd, ind in zip(rows, cols, inds):\n",
    "    InvM[rInd, cInd] = symbols('ogxm_'+str(ind))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Display the matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABFcAAAF3CAMAAABXFFReAAAANlBMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAABHL6OuAAAAEXRSTlMAMquZdlQQ\nQN0iRLtmie/NfG6M4XYAAAAJcEhZcwAADsQAAA7EAZUrDhsAACAASURBVHgB7V2LVuu6DgxQOFAo\n+/L/P3vtOLLs2HFla0qToq51jp2HxqNRKtrAzkxPP/PrebKXKWAKmAI6Bd5DO5mmp5+Xk3u96uAs\n2hQwBUyB6cs3k+cf31eeTA5TwBQwBWAKfFlfgWlpQKaAKRAUuNZXnj+hSr29XE4vyk9HIxjgNKqa\njPCqAoWdMDgIEAIEgZEIBoZj5CMAqzhqgoWxV/vKF+sNmP1zN3He3j9USCMYz9g0qvxHeFWBwk4Y\nHAQIAYLASAQDwzHyEYBVHDXBwthrfYXVRsy+vj3K5Z8GC4GhWX8rFswLBgcBQoAgMBLxwXCMfARg\nFUdNsDT2d/vK5/y77NefN65i9wyB0b2oIADMCwYHAUKAIDCSOoDhGPkIwCqOmmBp7O/2lZ+5r5x/\nNL/SRmDwRYSbgXnB4CBACBAERlIuMBwjHwFYxVETLI1t9pW3y8sp6P16cX/for6F+/bz4uHOPwsq\n11I+G8CopPHmfsP+cT5dHJHz5fPjyW1qPkM5+gO8WknD4CBACBAERqIYGI6RjwCs4qgJFsc2+8pl\nOs83RKaLaykfvim4d+JlLsHTSJP5+JmDn8LApeyZDWCUabw5Hl/vp+nd3em5TM/v7vOTT3GirGjs\nIDbAq4UOg4MAIUAQGIliYDhGPgKwiqMmWBzb6isfr9PLfIs13BBx316+3Pvx1X2XeXp5fudKiGdi\nVg3EfowyjeniPpxc3F2e14/JHf3nv505ISgrGhssykP9vEqMZA8MDgKEAEFg3EKgBHOegnkyPBBY\nBaUJFse2+op784VvLN/zz/Kfj+nd/+nJfNf1daSvvIUPKsrvQfNnng6MMo3Jf+f5Fz5x+aP+do9v\nMxNlRSNfFFdniNySRWBwECAECALjFgIlmPMUzJPhgcAqKE2wOLbVV9z7bO4hTz/+rz9O7+4egv/L\nk2//Nhx457mocH/lVXffdr5H04WxSsMR8VzoJs95TvLTf+OjrGicTxT+D5FbshQMDgKEAEFg3EKg\nBHOegnkyPBBYBaUJlsa2+4r7huB+in/N7eTzxf1DIv+D/t13mZF3ngudPyN8zTCsd99sAGOVxrxe\n+EdRPp+L/+gVbkhRVjT2MBvg1YKHwUGAECAIjEQxMBwjHwFYxVETLI1t9pU398HC31Fx/4Zo/raw\n9BX/g37knec+88y3gcNNGy5k36wfY52GY/82nXxSJ99X/vkPQF/ffkpZ0eh2iV/9vJrQMDgIEAIE\ngZFIBoZj5CMAqzhqgqWxzb5y/pk+3IeT+evPs/vQ8hS+B433lcn/5uXtW/cPhLox1mm4X3S/Tp+u\nr/jkXMN0H6KeAifqJzTyxSaYdfNqY8LgIEAIEARGohkYjpGPAKziqAkWxjb7yvR5mm9CnJ9PF/9t\n4W1+oMLP2ZVg6J3nENzfwTzr2soAxiqN6e35dHo7Xb7mtnL+Obt/Cxn+yRJlRSNfa4IZIrdkGRgc\nBAgBgsC4hUAJ5jwF82R4ILAKShMsjG33FZZkcrdX3P0R31LmuyND77wE7l7TOY1s8fn2yrKHsqIx\nO9E2TAFTQKiApK/4zxdv86eUk2su/u9X3P/DH8wJV9nFaZxGRme+vbLsoaxozE60DVPAFBAqIOgr\n8192hH9vNF1eLi/uBufHy/vPM/2eVrjSvU9L00i4+Ju6y4uyopH222gKmAJ9Cgj6yuvz1+klvvn6\n4Hd0dj2N1/effwdrkDvS1KiYAnUFBH2lHmh7TQFTwBTYUMD6yoYwttsUMAWGFbC+MiydBZoCpsCG\nAtZXNoSx3aaAKTCsgPWVYeks0BQwBTYUsL6yIYztNgVMgWEFrK8MS2eBpoApsKFAR18ROhJVFsK6\ngo3zWFMbRsImtKYVtofJreEgQAgQBEaSHBiOkY8ArOKoCRbGdvQVoSMR1yfOsK5g4zwioWUyjIRN\naE0rbA+TW8NBgBAgCIwkOTAcIx8BWMVREyyMlfcVqSMR1+c2MxwPHNINMoWRgwAhQBAYidBgOEY+\nArCKoyZYGivvK1JHIi7QbWY4HjikG2QKIwcBQoAgMBKhwXCMfARgFUdNsDRW3lekjkRcoNvMcDxw\nSDfIFEYOAoQAQWAkQoPhGPkIwCqOmmBprLiviB2JuEBhVnEFW5/Ssz3Mo1hkFKlMCO9ytjxvV2fh\nNmc8mmUmFwIEgZGQAsMx8hGAVRw1weJYcV8RO4dwhcKsdAUjc7PLy6d/6ELfa5hHscwoUpHQFZez\n+6Y5mmUmFwIEgZGQAsMx8hGAVRw1weLYm/eV0hWMzM0u7tmPz7PvGRf1+kyc2VWoQaQyobbL2Z3T\nHMwyVw8BgsBIWIHhGPkIwCqOmmBxrO8rb+7x9FdfYkeiFZKD3jA3e3d9JTyJexXS3BzlUYIOIpUJ\nefW2Xc7unOZglrlcCBAERsIKDMfIRwBWcdQEi2PdT1/35pY8yVrqSMQVWmYrV7BobvbtHpb7MT/e\nsghp7RjmUYCOIq0SmnG3Xc7uneZolplcCBAERkIKDMfIRwBWcdQES2PF34PGPcVWrmCJuZl7TG6/\nRZnUGYkvlK3ZKNIqIQ/fdDm7b5qjWWaqIUAQGAkpMBwjHwFYxVETLI2V9xWpIxEXKMzWrmCJuZl7\nwH//MyBHeax5DbukrRO66nJ23zQheiFAEBhJEcFwjHwEYBVHTbA0Vt5XRj3F1q5gibmZeww3l1M8\nEzojCfDGkNYJXXU5u3OaY1mu5EOAIDASWmA4Rj4CsIqjJlgY29FXhI5EXJ9ltnIFS8zNvkbayoAv\nWUFp2TGY0Tqhay5nd05zMMtcNAQIAiNhBYZj5CMAqzhqgoWxHX2FhR+e5eZmZ+f9fJbcMh5e7+aB\nhc1ZxeXsAdK8uY62wIMp8Gt9hV3ByNzs6d/r6+un4Ffc+5ScE0r5lS5nB08zTc7mpoBQgd/qK6kr\n2GJu9v3jX0KeuzstTYjJVVzOjp0mp2YzU0CuwG/1lbormJzn7s6sJmQuZ7urkxG6iwK/1Vfukpwt\nagqYAndRwPrKXWS3RU2Bh1bA+spDl9eSMwXuooDvK//9/HeXtW1RU8AUeEwF/uf6iuzfHT5m/paV\nKWAK4BWw70F4TQ3RFPjrClhf+etXgOVvCuAV6OgrQkeiJsdxjFtYgY2zqSU5jHaL1BKCw7wSjAkB\ngsBIOIHhGPkIwCqOmmBhbEdfEToScX0qs3GMW1iBjbOppDYNo90itYTgMK8EYzy5BARChPHAcMcC\nViWvCRbGyvuK1JGI61POEBgl6ugeLBss2mhOZRyEFwIEgZGkB4Zj5CMAqzhqgqWx8r4idSTiApUz\nBEaJOroHywaLNppTGQfhhQBBYCTpgeEY+QjAKo6aYGmsvK9IHYm4QOUMgVGiju7BssGijeZUxkF4\nIUAQGEl6YDhGPgKwiqMmWBor7itiRyIuUDEbxiitwArs/h3DbKpLjaJVUoPanI3yypJEgCAwElJg\nOEY+ArCKoyZYHCvuK2LnEK5QMRvGKKzApuhtdro8Dz7CZZhNkZffMYpWpnbN5qwv41FeWZIIEARG\nQgoMx8hHAFZx1ASLYw/RV0orMPI2e/mappdPvih6ZmKNRKCDaGVqU9vmrDfjQV55zggQBEbCCgzH\nyEcAVnHUBItjfV95Owl+4osdibhCxWwUw9Hb8DZ7cY+y9CmMvEbZ1NcaRCtTm3wxtm3OejMe5JUn\niQBBYCSswHCMfARgFUdNsDj2wxltyP59kNSRiCtUzoYxVlZg0dvML/HS7cW6EBtmUybm9oyirVIL\n2Ns2Z/54V8ajvAKR5f8IEARGQgoMx8hHAFZx1ARLY8Xfg8Z9ybhi4xgrK7DU2+zte/TR21KPpSSB\nxnQUbZXavELb5qwv41FeWaoIEARGQgoMx8hHAFZx1ARLY+V9RepIxAUqZ6MYaysw9jZ7O/3rtzZb\niI2yKfPyewbR1qk5c8S36eS/2s3fTufncH99++9Gr+9+GfettS/jQV7zWvF/CBAERiQ0rHeCsDEF\n8+RVgMAqKE2wNFbeV0Z9yVhXNxO6GmUxbmNtBZZ6mw3ftx1ms2a3bI/ltk7N5fo6fbq+8uFuSbtv\nV+6m9FP4RLb0Fbez7071GK8lKRoQIAgM4uNGMBwjHwFYxVETLIzt6CtCRyKuT2U2irG2Apud6hfP\n+Cf3Rhx7jbKprzaItkptertmc+ZviPVkPMgrTxIBgsBIWIHhGPkIwCqOmmBhbEdfYeHvOku9zeY7\nDR8/80/2u5LCLF64nE2FzdmDZYzRzVB2p8Ch+gpbgZG32fuH+z3zfPdhd8r2EeLUsrjS5uxhMs7y\ntI0HU+BIfSW1Alu8zZ4u7q9PXW85+itNLcmlYnP2KBknWdr08RQ4Ul+pWoE9RknqqZnN2WNU9w9m\ncaS+8gfLYymbAodUwPrKIctmpE2BXStgfWXX5TFypsAhFfB9xXzJDlk6I20K7FYB8yXbbWmMmClw\nWAXse9BhS2fETYHdKmB9ZbelMWKmwGEV+OW+InQ1asqJwGguYAdLBRCiIzASZmA4Rj4CsIqjJlgY\n+8t9RehqxDWuzBAYFVjb1VIAIToCI+EIhmPkIwCrOGqChbG/21ekrkZc43KGwChRbU9TAYToCIyE\nJBiOkY8ArOKoCZbG/m5fkboacZHLGQKjRLU9TQUQoiMwEpJgOEY+ArCKoyZYGvu7fUXqasRFLmcI\njBLV9jQVQIiOwEhIguEY+QjAKo6aYGnsVl/JvLGePz5Op+fTdL6cPs9cgO6Z2NWogYzAmKbL5XQa\nfSxuQu78/n2entzjX86Xz48np5F/XuTYi6DO7t9nK9XOaqflNWeDEB2BkUgLhmPkIwCrOGqCxbEb\nfSX3xvr8dM+Qffs5uZ7y+u0rQK5g01OXd4/YfYSrXMwgGP5hev/8Y/yJP43Fcs0db+5Zby8vz+6c\ny/TsMS9ejgXr8vL5Im8zEeoyFWpHRJkjWV67yCvyGcgVIvqP82Txz7qbBz9VvRCUqgSOAKziqAkW\nx270lcQb6+l1evdvwQDpz5/IFezp5Tk8y7laoHKnmFUZGvcgML79pf38OhF/GuMiwolvGx++rTh3\nMfdY/dnKiLAu7qkwz3IHEoKqqE2IUkeypHYJL+JDaMIUw2kI0REYCWkwHCMfAVjFURMsjvV9ouJL\n5i/z6I31Nj9INvjcPPufyO/+K8SPP4ef5ew2rr7ErkYNJADGy0x9XoP409hYuH7oybcTbyU2P3J2\nfjxTwPLPdQuP964HlntnqDf3udB/1UzUJpWljmRZ7SIv5jOQK0B0l9f8QeUcrZFKBXr2gOF46SMA\nqzhqgsWxDV+yeAG8+s8o08v80eR7/kbkn9D27Z/d3HmNSl2N/HpbLz3GN395I/40bi26tf/8shw5\nz73q039LDFjuxov7iNdxN4qg1mpHRL+S1JEs1i7yYj4juepFHzdu83lXXghKFVg4T14DyFgFpQmW\nxm58D3JisDdWuJjf/Y+b88/H9Ja4gnVeo1JXI65FOVNjLD13Rib+NJbLNfe8UlsJD7gOd7UYK3zo\naCLEgxFqrbY7IyJKHcm4dvODt5e7bcuHoIgW174+UYvulkBgJEzBcIx8BGAVR02wNHajr2TeWPPt\niPDx3F30H6/sCsZXPJelNZO6Gt0WY3l+f/p5a+S95t7vc1uZvwmV/mEuh3e5ZRpDrdV2OAs7oSNZ\nVruEF/EZyXUfhcsuCwSlDJA2jgCs4qgJlsbW+0rmjRU+zfsT3T2Xi/vVR+IK1nuNCl2NqMTVUY0x\n/wbn7cXfJCL+NFYX3Nr59P51+Tr9m7/r1PzDLvHTzBZC3M9QhdruHGYncCTLapf6mk0LH0aLy1+f\nqEV3SyAwEqZgOEY+ArCKoyZYGFvvK7k31nxj5TT/WH59vjy5O3Dzfdv59mLX74PcDU73dyPPyj8c\n0WNcLl9f4aMEvcdo5ItLMHNt4+s7tJXzj/vbnpdgDEBYX/K2MjHUea22I0KI/sup/5TVfOW1S3gR\nH0Zr4uQH9aJjip+wQlBK4Hh6BGAVR02wMLbeV1jj6uzdt5SB3wdVwe66k95jNI6SKfzD/J/Kuf+U\nLXShM7MbciRLeEU+2lxHNbK4P6TAUF8hVzD3k3T+M7nj6kX8aRzNpPQPe/r3+vr66X/nq38FdiOO\nZMyL+Whz1edjCA+vwFBfcd/UL/5vST9e3n/cn/cf9kX8aRxOpOIf9v3jX8OISSCxG3AkS3gRH0JL\nFrCpKYBWYKyvoFkcHG+v/mF75XXwchv96wpYX7mukZ1hCpgCfQpYX+nTy842BUyB6wpYX7mukZ1h\nCpgCfQr4vmK+ZH2a2dmmgCnQVsB8ydr62FFTwBToV8C+B/VrZhGmgCnQVsD6SlsfO2oKmAL9CnT0\nFaEjUT+Hzoi98OikfbfTd6MXmAgYjutzBGAVR02wMLajrwgdibg+N5rthceN0oPD7kYvMBEwHOt+\nBGAVR02wMFbeV6SORFyf28z2wuM22eFRd6MXmAgYjoU/ArCKoyZYGivvK1JHIi7QbWZ74XGb7PCo\nu9ELTAQMx8IfAVjFURMsjZX3FakjERfoNrO98LhNdnjU3egFJgKGY+GPAKziqAmWxor7itiRiAuU\n+WOp3bYW3AEexAho/UWQYQS5nE1MEOImNpMb1SsrHoLOKJFc6rgFhou47qll8/O4UL4BNwFWcdQE\ni2PFfUXsHBKFzP2xCrct8jaLflkxsDnp50FwEusvsgCjGMn44R/Nl7qcdWYU10gIRjcxInSReZJF\nrDgZ1CsvXqRDqdEYl7k+GSSyBQyG42WOAKziqAkWx96wryT+WBW3LfI2I78sLmx7Js6sgPHPWPIu\nYhUyZNZFYxHb2rF2OevNKGITwcRNjAhJPckiVpwM6pUUL6FDqdEYVxFMBolsIYPheJkjAKs4aoLF\nsb6vVHzJWGaaiR2JKMCbdQm8zdgvKwY2J/08Erhr1l/po2STsOa0cDnrzShFXwhGNzEiJPUkS7HC\nfFCvrHiRDqVGY7nc9p5BIluAYDhe5gjAKo6aYHFsw5eMpZ5nUkeiLCz6Y63dtt6cEdHsbcZ+WVng\n9sYQjwAnsP7qf/hr4XLWnRHnSgSjmxj1FX+K1JOM4ebZuF6xeJEOpUbjaqn25jiRKi4Yjtc4ArCK\noyZYGiv+HjRkK8X+WGu3rcTbzL132NiUy7s1kzojlfES66/uvrJ08Hk1Du7KKDKNBBM3sYgp9SSL\naMtkWC8uXkKHi9Wd4jCRdUZhGwzHixwBWMVREyyNlfcVqSNRLFDmj7V220q8zcgvKwa2J908CE5i\n/ZVYalDYtbHictaZEa3ABBM3saWvCD3JCCoZx/TKipfQ4dQ6TNcCmzEiSSb5FAzH4EcAVnHUBEtj\n5X2l11Yq88cq3LYSbzPyy+LKtmdCZ6Q1iMz6K348WIdvbldcznozCthMMHUTY0ICT7IqyRG9suKl\ndGJqHaZrxGqECMVWRjAcr3AEYBVHTbAwtqOvCB2JqD65P9babSvxNiO/LAq8NnbyIDiZ9Re/jSnu\n6li6nPVmFJZICCZuYkxI4ElW5TqiV168hA6lRmN1xY2dI0Q2oPxuMByvdARgFUdNsDC2o6+w8IhZ\n9DaLflkI1GEMevvSOAS0BKszStzEwhezIU+yoRRqQQkdSo3G2um2zxSY7tZXyNuM/bLuWg0y66Jx\niEwI1mfEbmLuNuns/OZ/sfv17X/5e4cX06HUaLwDGVvyEArcra+Qtxn5Zd1VLTLronGIDAWrM0rc\nxAhzwJNsKIdaUEKHUqOxdrrtMwXcD0H3d3HhF4qmxk4U2Jmb2M7o7KRIRqOtgPWVtj521BQwBfoV\nsL7Sr5lFmAKmQFsB6yttfeyoKWAK9Cvg+4r5kvXrZhGmgCmwrYD5km1rY0dMAVNgTAH7HjSmm0WZ\nAqbAtgLWV7a1sSOmgCkwpkBHXxE6EjV5IDCaC9jBQgGI5hCQgppqx80oHQFYxVETLIzt6CtCR6Lm\ntYLAaC5gBwsFIJpDQApqqh03o3QEYBVHTbAwVt5XpI5ErWsFgdHCt2OlAhDNISAlN82em1E6ArCK\noyZYGivvK1JHota1gsBo4duxUgGI5hCQkptmz80oHQFYxVETLI2V9xWpI1HrWkFgtPDtWKkARHMI\nSMlNs+dmlI4ArOKoCZbGivuK2JGoca2MYGT+WHpzM5z1FyFN54ueViYayuZsBh3RPGPjNwZBstoh\nvM0SYoOUEoSN6RGAVRw1weJYcV8RO4dslMvvHsDI/bG2zc2Etl0t6y+yAGskkByKSNNlKmgR1IB9\n11TanAlzS8gl0wHNk+hlOgaS1670Nrt/XmWmQxdpDabcNyZiieP2qKA0weLYvfeVxB+r4idG5mZi\n266G9RdZgFULWe4kpJbN2Yh917S2ORPnVnL0e8RXQj087B0DSWpX8TbbQ161nMdyrSGt9gGBVVCa\nYHGs7ys38iVbyeo2xa5GHJr5Y739nN2RYDDx/Omm70/uf84kpMe2a8v6ywHPD+Hlta/NZiRnv1bQ\nIqgR+67C5qwntwrjAc1LlDGQrHaFt9ke8iozHbpIazDlvjERSxy3RwWlCRbH3tiXbKWK1NUoD4v+\nWJvmZv58oW1Xw/qrs68Qkmt08/pzV/o+zdwD1Ih9V2Fz1pHbvPT6f2Oar1CGQWLtCm+zXeS1SnPe\nHM61BpbuAwKroDTB0ljx96AhX7JUVT+XuhplceyP1TI3E9p2tay/+vpKRFo62vvFsT47G0f/Y5qh\nOu27lh8IswARRJhbplrcGNI8Ri+TURCuXc3b7P55rfP026O51rCyfUBgFZQmWBor7ytSR6JMytVG\nP0bmj7Vtbia17WpZfyXNYEW7tslI4Y5I+DLkOp+7lZD2lV77rtLmTJpbjaXb1695BWgIJKtd6W22\nh7wqqWIEuzHwUD2IkyZYGivvK72+ZJRFNgpdjWJM5o/VNDebJLZdbeuv+PkgLr89SZAKWi4qQnXb\nd9VszkS5bXLt1bwKNACS1a7qbbaDvGrJDuRagyn3AYFVUJpgYWxHXxE6EpVyJnt6MXJ/rIa5mX/6\nt/+g0H61rb9iM2iDzEcZaTqvabkTCGrAvqu0OZPltkm6V/Mq0ABIXruKt9ke8qolO5BrDabcBwRW\nQWmChbEdfaXU6e57FnOzAduuxGsrNgFqBoC0FiidfdcMMpAbgD8cItE7iPIgecGFehDAY/cVMjfr\nt+1iry334WK2/uIRUFqIQ9kOLMkAUgQI1ps8zfprBiNjQDdX4Nh9hczNum27Eq8tsv6iESA5Qans\nuwikOzdAAnCIRG8S5SHyggv1KIAH7yuDZTCvrUHhBsNM70HhDhv2N/vKYctlxE2BQyhgfeUQZTKS\npsChFLC+cqhyGVlT4BAK+L5ivmSHKJWRNAUOo4D5kh2mVEbUFDiMAvY96DClMqKmwGEU6OgrQueQ\nZup7wWiSHDyIyC1ZGgyXINu0rcDNlAcCq6A0wcLYjr4idA5p1mwvGE2SgwcRuSVLg+ESZJu2FbiZ\n8kBgFZQmWBgr7ytS55BWzfaC0eI4egyRW7I2GC5BtmlbgZspDwRWQWmCpbHyviJ1DmkVbS8YLY6j\nxxC5JWuD4RJkm7YVuJnyQGAVlCZYGivvK1LnkFbR9oLR4jh6DJFbsjYYLkG2aVuBmykPBFZBaYKl\nseK+InYOaRTtrhiZmQ3Y8cfljMgtkQ4CpzIiSvyRxqx/EgCM3FkFxzglCm9NIcrXwIHAKihNsDhW\n3FfET/ivSbrsuydGbmZTOv40aIsOIXJLFgLAsRFRgiuepv5I0fpHHO1OTAEgcucVHOIk4Q9Qvr4M\nEFgFpQkWx/6ZvpKY2VSMiNyjy7xviHvRGLbk/xcrLoMEwJERkXuYxOeLf5h3V24+4OPZ/+91+udH\nf6UQEo1u9+YrAtTkJoBLhzdZUsGME0GMVm6VAUD5FeKyCQRWQWmCxbHiviJ2DqmLOu+9J4a/zP+F\n3lFx/CFTMhobSWwcQuSWQOvhohERuaN150b+SNH6h5BoTPhWpwRQGCwRQJc3WVbByIkgurOrMnY7\n9cpvIAOBVVCaYHGs7yuyfx8kdQ7ZUHXefWeMaGaz6fjDT6ZtpVE9hsgtAVbDRSMidkfre9Bm9EeK\n1j+ERGNCtzaNAIXcBNDtTRYrGDkxRF92NcJhn1r5LWggsApKEyyNlf/7IKlzyJaqfv99MdjMZm1E\n5KjRVUljK43qMURuCbAWbvnJ4hDZHa0rN/ZHitY/hERjQrcyZYBC7hQgHKzEV3ZxBSMnf1aA6Mqu\nAk67tMoTTjECgVVQmmBprPh7EMRWReo+UpQk2TGIkZnZrI2IHDxdlTQmK8qmg7y2wNVwiRGRy86Z\nzXKOW2um+xN/pMz6h3zWaExjsnkCUJGbKE1yb7KsgimnBWK4chltkOHSCnPeVJeUQVVQmmBprLyv\n3MU/iJXkmdDBhAP8LDOzaTn+jF+dQ7xylumWFi4xInIfFGd7147cEn+kzPpnQSLElHA+TwBqci8A\nHd5kWQUTThGiI7uc6npLq/waL24DgVVQmmBhbEdfETqHRBVrk7th5GY2Dcef8asTkVuimRqOjYjc\n73Fm4I7cUn+kxPqHkOKYEM6nKUBFbgaQ+Ml56LyCKafF3qwju5zqekut/BqQtoHAKihNsDC2o6+Q\nOo850lVJ4yNlSe5oY7kl1j8TIdE4rBIDSPzkimVSTmRvNpZdAW07EApYX1lUpKuSRoS4O8GI7mhj\nubH1z0RINA4nqPUmi5wSe7Ox7IZTsMCWAtZXFnVuYE7W0v0Xj5ERmLtVuhiwdS2eWP8QEo1dOOnJ\nBOB/3/z17W8p970STgwxlF3funa2VAHrK7NSZAJGo1S+I5xHRmBjuaXWP4RE43D2BDDoTZZyIoix\n7IYzsMC2AtZX2vrYUVPAFOhXwPpKv2YWYQqYAm0FrK+09bGjpoAp0K+A7ytvp/47Z/0rWYQpYAr8\nFQU+3B9ihn918VcytjxNAVPg1grY96BbK2z4BlfIfgAAIABJREFUpsDfU8D6yt+ruWVsCtxagY6+\nInQkajLeC0YgiWDTTFdzcNfkNImhYo8nEJCxCkoTLIzt6CtCR6LmZbMXjEASwaaZrubgrslpEkPF\nHk8gIGMVlCZYGCvvK1JHotZlsxeMwBHBppWt6tiuyakyAwUfTyAgYxWUJlgaK+8rUkei1mWzF4zA\nEcGmla3q2K7JqTIDBR9PICBjFZQmWBor7ytSR6LWZbMXjMARwaaVrerYrsmpMgMFH08gIGMVlCZY\nGivuK2JHosZlsxeMQFHHJjPJwvhupcqNkwPbgam8zSYic76gfcTGBUplXq6E0+n543y6uL/lwhON\nywEZq6A0weJYcV8RO4dEHcvJXjACMxWb3CSr9N1azGzIJ6fU4sqeYXJNO7BIR+q1E73NKJLGK/SX\nw5HMZYo+YoTQYx1UWW1YoBIrryUTXYyJiHAZ2LcHyFgFpQkWx1pfefq59F0f89mJSVbFd4vMbMgn\np38BcQHX0P6fZMx+YhVaRIforWPLbfI2o0gayzOre4hM4iNGCOT7U40T7BwWqMROapkQJYJEuIzr\n3ANkrILSBItjxX1F7EjUUHsvGIGiio1/x2zbnNGj78knp6HJxiENuU07MKYjfLZa9DajSBo3WFd2\nz2ScWj+v7uDl540Q2PenEiTYpRFoBZ/VMhIlgkR4FdS/iWQcfhaeo51SDxsND3Gs7yt/xpdspb7U\nY2kVFjdjVQvfLeorqU9ODJNNxslt24ExHWFfid5mFEmjLAd3FpGJPmIpQo91ULHiuEAFlNsRaxmJ\n+rM8wZRwLVK+D8hYBaUJlsb+JV+y9RUg9Vhaxy3bbJIV3h/v/vvU+edj8j//+I171WinDj9MrmEH\n5lfqcRJafjotBCkRGuu8872RTOYjtiDIrYNy0LA1LFANjGuZEo0Ee1KuwcMZq5LXBEtjxd+Dju5L\nVim31GOpEjplJlkV3y3uK8G6p4bR3DdK7podWJeTUOZtRonQ2KQfDjKZ1EcseBBF3x8BTu2UUYFK\nrKyWTDQh2JFyCR/34Bjr3owaHtJYeV85ti9ZLG46EXospSHLPDPJqvluxb6yWPdUMK7sGiN3zQ6s\nz0ko9TajRGi8Qt8fTsj8fLrN76c5KCJIrYPqS40JVGJltUwMz9yZC8FIuAzu2oNi7BZVQWmChbEd\nfUXoSNRUei8YgeQ4m9wkq+K7RX2FfXKaulQOjpG7YgdGdIheZd1sF3ubUSSN2WkbG0zmnPiIMcKQ\ndVBca0ygGB4neS0TomRMxIRjzNgExditroLSBAtjO/rKmJZ/NGp546qNdrDyRTrSvhKXp0ga4wHZ\nJPERCwiJ748M4ZfOikSZ4GDKv0R4p8tYX7lNYYKZDfnk3GaNblSm0+u1Q5E09i4dfcQmQvC/vh2x\nDupdufN8JkoEiXAn0B8/3frKLS4AMrMhn5xbrDGASXSInhyCImmUR85nJj5ihEC+P51INz49IUoE\nifCNV34weOsrD1bQPaaT+ojtkV/kdBiikfFeJ9ZX9loZ42UKHFcB6yvHrZ0xNwX2qoD1lb1WxniZ\nAsdVwPcV8yU7bv2MuSmwRwXMl2yPVTFOpsCxFbDvQceun7E3BfaogPWVPVbFOJkCx1ago68IHYma\neuwFI5BEsOF0sWiMa7O6AjfT+2bA9TyG9qo4aoKFsR19RehI1JRpLxiBJIINp4tFY1yb1RW4md43\nA67nMbRXxVETLIyV9xWpI1FLpr1gBI4INpwtFo1xbVZX4GZ63wy4nsfQXhVHTbA0Vt5XpI5ELZ32\nghE4Ithwtlg0xrVZXYGb6X0z4HoeQ3tVHDXB0lh5X5E6ErV02gtG4Ihgw9li0RjXZnUFbqb3zYDr\neQztVXHUBEtjxX1F7EjU0GkvGIEigg0nq0MDu5yRHdiwzRYBnC8QyzWduRmLnMx0eidA6ykSOCsr\n0JlNxVETLI4V9xWxc8i6Tsn2XjACJQQbTk6FljtjbbucLU5ZvGp9Fu3Apoof2CTwJIsAl6kgQzZd\nNNYppHtLczNhHinIeq7Sew2WbgOB87KWhmfDMqg4aoLFsdZXxnzJ0itxnosVLyLdjsQZq2InRjZi\n5JRVQ8j2eUOA2Zsssdkiey0CywLWGwRQIUM4NK5DK9trczNxHhWsuEuld0SpTIDASVmTSlD6NFY4\nXNul4qgJFseK+4rYkaihyV4wAkUEG05WhZY5Y739nB1ssJZ4do+b9hvzE3TJKWvedeV/5E0WbbbY\nXkv2DMrFT6wgQzg0XuHhDhfmZj15bMKr9N5EdQeAwFlZYyUofRpbZDaOqThqgsWxvq+YL9lG+bp2\nSx2btkCjM9amy5mPFFp5kR2Y8zPyF/fnd2KvJeorBFCQIZsuGrfS4f2FuZk/JMyDUYqZVu8CkHZg\ngWNZYyX8OpQ+jbS2dFRx1ARLY82X7Gt+50kLun2e1LFpA4GdscK1Vnc5i05ZGyjL7mgHNqU2Wx2e\nZBGgQoY+S/HYJLP8iCNmQW1hHi1gpd7b0FBgLmtaCUqfxm0yG0dUHDXB0ljx9yCdFdIij9TVaEPN\neTcCI+DjkDyeBi1zxtp2OUucskIGG/9nO7CJbbbcuWJPMgaokCEcHjdoLLtLczNpHk1cjd6/BJyV\nlStB6dPYZFM/qEpeEyyNlfcVnRXSIo/Q1aguJhAjQCHYMNVxtMwZq+lyRk5ZvGplltiBZTZbi73W\n9e9BDFAjM5FNF40VCumumrmZKI8UpDIf17sClu5CAWdlzSoR0x92ZlNx1AQLYzv6itCRKK1QMd8L\nRiCGYMMpjqPlzlgNlzNyyuJFazO2A3O3V86X04uz1HAvste63lcY4FwhQzg01ihk+0pzM1keGUi5\nMa53iZXtQQHnZU0qwekPO7OpOGqChbEdfSXT3jZ+U4G5Fcxfxj+WbxXC1aPNlnOkd7b0Z29ker2v\nNMEJh8bmyfnBEDKSR45zxK1YCUqfxiMmc52z9ZXrGt3/jGAjRk5ZHXzYZovttXo9yfLlCIfG/Ghz\ni0IG8mjiHuIgV4LSp/EQ9HtJWl/pVez3zycbMXLKkjNIbLbIXovA5CD5mYRDY360uUUh/Xk0YQ9x\nMKkEpU/jIfj3krS+0qvYgc43m629FOvPVcL6yl4uPeNhCjyOAtZXHqeWlokpsBcFrK/spRLGwxR4\nHAV8XzFfssepp2ViCuxBAfMl20MVjIMp8FgK2Pegx6qnZWMK7EEB6yt7qIJxMAUeS4GOviJ0JGrq\nsxeMQBLBhtPFok1gOOY5NEOwQWAMkbegXAFNIYSxHX1F6EiUp7Da2gtGoIVgwwli0SYwHPMcmiHY\nIDCGyFtQroCmEMJYeV+ROhLlKeRbe8EIrBBsOD8s2gSGY55DMwQbBMYQeQvKFdAUQhor7ytSR6I8\nh3xrLxiBFYIN54dFm8BwzHNohmCDwBgib0G5AppCSGPlfUXqSJTnkG/tBSOwQrDh/LBoExiOeQ7N\nEGwQGEPkLShXQFMIaay4r4gdifIcsq29YARSCDacHhbNPRXePVzJPTAlPnaZV2rOMicsiKeYX26U\nTUoVgbHg4W3O2IcN6B5GoO7RN/pKZIVVkdQUQhwr7iti55D0WlrN94IRaCHYcIJYtGkQLnfC2vQU\nk3iTcWrOjOjHPRTKP+VsHtIj8jkCY16Nbc7I10vgtNbmmfiwle5hL58v3tOg+xVBnTtcWYmF+4Vy\nuAKfFzaSdFHBCaZDAk0hxLHWV1RvFr4axIpzSGs2CJc4YTU8xUTeZCm7QTYpBKQ3zYBkc0a+Xt3Z\nZLTmDd83vJFbxT2sw39thUugU6USxJ3GVWi5mRQ2IenOm32muiTQFFMcK+4rYkeiUpO4Zy8YgRCC\nTUwNaWflQQfJ+Wv5X/gJ9uYwVgZn/gllTz/+abedz6IcZOMTiS8EhgeLNmfs69WZTaSUThYftsI9\njDVLzxbOyR2uqARxp/EqXlbYSNJdJ5dgnNIhgaYQ4ljfV8yX7GpZBSdIHZsEUP6Ucbh4S2bTU6y7\nryjYcLrjGTGGm0WbM7832Bt1vKkyqGSDfNgK9zC5/1qCtkwJ1H2o8Hte5oeQf5+Wo+RJRuOye3uI\nhU1Int66+4qqmNIimi/ZTnzJ1peT1ABqHec+kPiHY/ufbw1Psd534jCbhB4CI36QC7jzs6f7u2RC\naplGH7aae9ji51ZGXdkTQauVWLhPNF4B8580qbBM8ul1oK9oCiGNFX8PUvlukWZSVyM6vzYiMAIu\nDsnjYdEG4TInrIanWG9fQSSHwHA6s81Z9PXqzSaUP/k/+7CV7mHutODnlpwvmjLoVFaCuNN4DTEr\nLJO8uC+7/kdI1xdbTSGksfK+Yr5kV2ovdGy6ghIPj8BlTlgtT7Hud+IIm5jKMkFguE9h7v6q+8dT\n/me325jvJXVnM8fy/9iHreoeJvRfY7x5loBWK9HlTZYVlkl+fYz0FdU7WVjEjr4idCRayZtv7gUj\nsEKw4fywaO6O3On0HN49vMaVWe6E1fAU634njrBZk0VgOEy2OSN/r+5sVswSH7aKe5jYfy1HZdCp\n5u5G3HnMw/OtvLBE8s3dqhn4vDJ0aREfYRE7+goh23hYBaKXmPaduAMFEl8vXDaFe1ji5wbMmbjT\n2AcdSb6+uNeP/1111/egvtXGzra+MqbbIaPIGMxdht+HTCAjzb5euGxK9zDWLFtcuUHcaeyCY5Iu\n7GO5v7Kzglpf6SrpsU8mYzCtN9k+VCBfL2A2Ffcw0gybM3GnsQc9IelMtz9/nl8noAQ9TBrnWl9p\niGOH/pYCh3APOwTJyfrK33rrWLamwG8oYH3lN1S2NUyBv6WA9ZW/VW/L1hT4DQV8XzFfst9Q2tYw\nBf6OAuZL9ndqbZmaAr+lgH0P+i2lbR1T4O8oYH3l79TaMjUFfkuBjr4idCRqMt8LRiCJYMPpYtFw\nvmRgXpxw72w3RHqJ7/D8u2kpXLijrwgdiZpF2AtGIIlgw+li0XC+ZGBenHDvbDdEeonv8Py7aSlc\nWN5XpI5ErSLsBSNwRLDhbLFoOF8yMC9OuHe2GyK9xHd4/t20lC4s7ytSR6JWFfaCETgi2HC2WDSc\nLxmYFyfcO9sNkV7iOzz/blpKF5b3FakjUasKe8EIHBFsOFssGs6XDMyLE+6d7YZIL/Ednn83LaUL\ni/uK2JGoUYW9YASKCDacLBbNPa5nzJeMCS2zEaDMA0tvqYWRm1y+ECZfhUrLjoMYnnVfHVlBNaZm\n4qtJ3FfEziFbNXP794IRKCLYcLJYNIhUM7kBXrkHVmmptRh1XfoMuwaIsLruj8I/T6cX/wzKhslX\nr91aiu/mheFZZ4YrOL+ZsI5eYmRFRmMl7NquTi3zgkYibpX5MZ4daYoXtr7yUL5k5QUpvhI4NPHA\nqlhqkVEXjRzXng0QSQH906G9dViFEZl7dblzpdjLfG141pthBZJYJ15ixJbGStTVXZ1aJgVNiLhV\nZlOznjTFC4v7itiRqCHKXjACRQQbThaL5j7pzpal3f7MTGiZDQD5N8N1c7New64BIqtkFuuwTZMv\n5dMYC8Oz3gxXdJfNLcMzsSVZBbZTy6yghalZT5rihX1fMV+ySum6d0kdm4TAMLgxoOiBVVhqkVEX\njcJ8NE5rYQly+SoY+cMAg7LC8Kw7w5oUxDrxEiO2PNYCm/v6ixoLmhAJpmZdaUoXNl+yR/MlW1+P\nUiepLI49sBrmZsvn6CywsTFEhPGiy1eN0WLupXmA9vKjOCy44HVmyGTjLLJmLzF31+U7+CzQGM8W\nT7q15IIykWhq1pGmdGHx96BBo6xcKamrUR6VbyEwAiIOyeNh0XBw/bwyD6zSUsvlSkZdNAY9r/y/\nn0gKyC5fJaNo7qXpKzXDM840pdIxZ9bsJUZsaeyAi6f2aZkVlIlEU7OONKULy/uKys2IBBG6GtHp\n1RGBEYBxSB4Pi4aD6+WVeWBVLbXIqIvGapnKnb1EUgR2+aoyQhiUVQzPps4MU8Z+zqzZS8zvD3Zq\nPPp9fa8eLbOCMhE2NetJU7hwR18ROhI19dkLRiCJYMPpYtFU5lFMys16eeUeWA1zs17Drl4iaRbs\n8tUy+VJ9XqkYnvVmmDL284Q1eYnNpzz9vGbjvNH1vx4t84ISkcTUrCdN4cIdfaUrbzv5YRUgczMa\n753ofJfiY7Zt1vWVJRHGA2YYvcQIncbfVi8SYVMzYJoxG+srUQqbiBQgoy4aRUE3PYnNvTAGZYSH\nzJC9xAidxptKU4IzEXdsNjVDphnXs74SpbCJSAEy6qJRFHTTk8jcC+XORXjADBMvMUKn8abSFOAJ\nETI1A6bJy1lfYS1sZgrcRoHdeIn9GhHrK7e5kgzVFPjLClhf+cvVt9xNgdsoYH3lNroaqinwlxXw\nfcV8yf7yFWC5mwJ4BcyXDK+pIZoCf10B+x70168Ay98UwCtgfQWvqSGaAn9dgY6+InQkaiq6F4xA\nEsGG08Wi7cyXDJwcq2azIQVU9VAFy+h29BWhI1Fz3b1gBJIINpwuFm1nvmTg5Fg1mw0poKqHKlhG\nV95XpI5ErXX3ghE4Ithwtli0nfmSgZNj1Ww2pICqHqpgIV15X5E6ErUW3gtG4Ihgw9li0XbmSwZO\njlWz2ZACqnqogoV05X1F6kjUWngvGIEjgg1ni0XbmS8ZODlWzWZDCqjqoQoW0hX3FbEjUWPhvWAE\nigg2nCwWrdt5iomsZhBeAyCZFRbK24y9yTTmWiuFGBRFc7WA2wQbng3Ugzl1BmeFFMsu7iti5xBO\noJjtBSMQQ7DhFLFod/Ul46Ro1p9cboW17W12ujx7FwrZq+XyNexJloCWNBd6Hc5dtUxKw7OutEvI\n/nokGH3BeSFLT7ONTKyvmC9Zcs1tTPsuRQ+SWGFVnMTICqvTm6vh8jXuSUagFZpEj+huqHN199rw\njHCvBm6d0F+PBKkvOClkxdNsMxNxXxE7EiUZrKd7wQi8EGw4QyzaXX3JOCma9SeXWWEVTmL+WWlP\nPx9TvzfXlsuXypPsmuEZ0SU5OsfC8Kw/7dWK/fVIAPqCs0IWnmabmfi+Yr5kierDU6ljk3ABGBwE\naAgkWmEVTmKpFVawAxKq0nD5Gn+2LYEWND0pTy+lKySanlYYnhFuelLffKgetER3cCxk4Wm2mYn5\nkpkvGV1vjVFqR5VCsBVWzUkseAOzR1cauTlvuXwN95UIWqM5P97a85mtjDeJtQ4snw/CKYRHYyuw\ncWykHhGuN5gLWfM0q2ci/h4EMcqSuhpFCSoTBEaAxSF5PCwaDg7Cqxsks8IqncScXrOnWZ83V8vl\naxrtKwxa0kzodVmwheuL/j9bBbjO5J9I8u/k99JIZ3SP3fVIV+gKzgpZepptZSLvKxCjLKGrUapC\nMUdgBFAcksfDouHgILw6QTIrrKqTWHT8Io+uos7FjrbL12BfYdAqTbIQi3QLVtd31AzPCPd6dP2M\nznrkIB3BWSGrnmb1TDr6itCRKE9htbUXjEALwYYTxKJ124kxkdUMwqsTJLfCanib+fu3waNrxbqy\n2Xb5GuwrCWiFJtHrce4qmV8uX1/z5xR3iNKlsTxbtKezHjlmR3BeyIqnGWeUrdHRV7I42zAFFAoE\nK6z5q3lwFOvDiuZaCcJgX2kszOAY5y7Co7Gx9D4PRdnJ02w7E+sr+yzhY7MiK6xRby4212IEjCdZ\npjuBE93s4MAG4dE4AHHXEJbd0Zg9zTYzsb5y10r90cXJCmvQmysx1yIElCdZVhACJ7rZwYENwqNx\nAOKeIYns5Gm2mYn1lXtWytYeUeDXzLVGyD1uTJfs1lce90KwzEyBeylgfeVeytu6psDjKmB95XFr\na5mZAvdSwPcV2b8PuhdDW9cUMAWOpoD83wcdLTPjawqYAvdSwL4H3Ut5W9cUeFwFrK88bm0tM1Pg\nXgp09BWEm9FeMILcCDZcOCzaA/qSgQVi5cGzm/G8GXCnABoewtiOvoJwM9oLRigEgg2XFIv2gL5k\nYIFYefDsZjxvBtwpgIaHMFbeVxBuRnvBCHVAsOGKYtEe0JcMLBArD57djOfNgDsF0PCQxsr7CsLN\naC8YoRAINlxSLNoD+pKBBWLlwbOb8bwZcKcAGh7SWHlfQbgZ7QUjFALBhkuKRXtAXzKwQKw8eHYz\nnjcD7hRAw0MaK+4rnW5G1VT3ghHIIdhwmli04/uSsTLLbEygy8vp5TwjvF5Op9dP/wzH0/PH+XRx\nD0o6X27gJDbGs0i33KEGLrUoFxHs0fAQx4r7Sp/rSD29vWAEdgg2nCcW7fC+ZCwMzUYE+nh3PeXN\nGYI4QyLXUj7ck+Rzn6xNJ7FhmzKc8pQ4jSMCUKwbK1r4o04W93oKwzy/9j8ND3Gs9RXzJbt2Ifo3\n9MWfpJJqBGN+kvX07hYPj8N3j6xMfLIaTmLjNmWQXKuKjgiQAFW0cEdnXfqy1fAQx/q+8nby5kNX\nXn1uRnWwvWAEdgg2nCcW7fC+ZCwMzQYEWgy9/Hf62YPHW5llPlmF4Rn7ZA0/lnKAJ6XYHnXANS3c\nW/fyM793e7LV8BDHfrivqcEfpC3KNHW7GVUA94IRqCHYcJJYNIjcMzkILwRIP0Yw9PI/I59mr4zT\n/GRrdynSQ6g3ncTcD/LlXC6QdNbPU4isAq5rcXrr7yuqS0uag/h70NTrZlQTey8YgRuCDWeJRYPI\nPZOD8EKAdGMsPxwvP0/T13yP5dM9PN+92Cer4SQ23le6ec6kBP/TANe1eHod6SsaHtJYeV/pcjPa\nUHkvGIEegg0nikU7ui8Z6xJn/QIFQy9/YyF8MPGOIJlPVsNJbLyv9POMKbYnKuCaFtPFfV/u/h6k\nurSkOcj7CsQoq8MRabNGCIwAjkPyeFg0HByEFwKkG2M29Jo/pMy/E3p2H1oyn6yWk9h4X8Epv76E\nuwVIACpaTF/udtNAX1ElKMyho690uBkleuTTvWAEVgg2nB8W7di+ZKxKMusX6HKhv145P58u3r4m\n98lqOIkp+ko/zyTJ1lQFXNHC3WYa6isaHsLYjr7SUsyOmQI3V2C5vdJYJ/HJUvSVxgK7OeS1IHcw\nR2p32Vpf2c2lYkS2FXhyh8IvlbfP8UfYJ+sGNmXtpX/taKbF7A7m+sr3ry0vW8j6ikwnO+ueCsx/\npRH+yVubBvlk3cSmrL30bx3NtPj6/Hl+nXaYrfWV37oebJ1xBV6fv04vUnv48WUOEXkMLayvHOJi\nMpKmwKEUsL5yqHIZWVPgEApYXzlEmYykKXAoBXxfMV+yQ5XMyJoCu1fAfMl2XyIjaAocTgH7HnS4\nkhlhU2D3Clhf2X2JjKApcDgFOvqK0JGoKcFeMAJJBBtOF4tmvmSsbH0G1ru+CHQvkDEQqi9F4cId\nfUXoSNSkuReMQBLBhtPFopkvGStbn4H1ri8C3QtkDITqS1G4sLyvSB2JWjT3ghE4Ithwtlg08yVj\nZeszsN71RaB7gYyBUH0pSheW9xWpI1GL514wAkcEG84Wi2a+ZKxsfQbWu74IdC+QMRCqL0XpwvK+\nInUkavHcC0bgiGDD2WLRzJeMla3PwHrXF4HuBTIGQvWlKF1Y3FfEjkQNnnvBCBQRbDhZLJp7KMD8\nLNdzfEQ0r9Q3gwAhQBAYSepauIrJ140dz7SMdclX8k0ApVNxDuK+InYOaVDcC0agiGDDyWLRcC42\nEF4IEAQGy60UqGbyJXQ8u7x8vgh8cRKqyxQoQDdULV9Ha3Yz60lIvLD1FZXZFl89YsU5pDWDwUGA\nECAIjEQxHVzN5EvmeHZx5ovP/xIi4qmOcbZMN1QtX/c0KP9w3K6ExAv7vmK+ZFnRxjbEjk0yeBgc\nBAgBgsBItFPBVU2+/IeQf8GP9G15OF1wWXz2O8nxzD+Rzvuj9b9UjPPleqGq+S6mZl0JiReefck8\n9NWX1JGoBbQXjMARwYazxaKpzKOYlJtBeCFAEBhJZhq4usmXF+uq49m3c4wOPgAJF9lUw3i1QidU\nPd9gataXkHRh/8Bdmd+h1JFopUC2uReMQArBhtPDopkvGStbnyn0Xn7qFoZnMsez5ftDnVVrr4Lx\nGrYPqp5vNDXrSUi6sPj+isrNiGSRuhrR+bURgRFwcUgeD4uGg4PwQoAgMJILQgNXNfmSOp65p3PT\np5qEjmCqYbyC74Sq5htNzXoSki4s7ysqNyOSRehqRKdXRwRGAMYheTwsGg4OwgsBgsBIrggFXM3k\nS+x4Nl2CnWtCRThVMF6v0AdVy5dNzboSEi7c0VeEjkRrBbLtvWAEUgg2nB4WzXzJWNn6TKN3afIl\ndzz7Gm0ruJK637VcTqdn7/ghe5X5vrGpWVdCwoU7+oosAzvLFDiaAtcNzyZ2PDs7u+iz/B29Ry1y\nU7ObJGR9ZY+FN06/pUBm8tValBzPnv69vr5+Dv1hXAv+t45l+c6mZrdJyPrKb1XU1tmhApnJV5Mf\nOZ59//hX89QdH8zyDaZmt0nI+sqOrwKjdmsFjmHyhVPh1/K1voIrmiGZAqZAUMD6il0JpoApgFbA\n+gpaUcMzBUwB6yt2DZgCpgBaAesraEUNzxQwBayv2DVgCpgCaAWsr6AVNTxTwBTo6CtCR6KmpnvB\nCCQRbDhdLJr5krGy9RlYb17kZsC8hHqm4qgKllHv6CtCR6LmunvBCCQRbDhdLJr5krGy9RlYb17k\nZsC8hHqm4qgKllGX9xWpI1Fr3b1gBI4INpwtFs18yVjZ+gysNy9yM2BeQj1TcVQFC6nL+4rUkai1\n8F4wAkcEG84Wi2a+ZKxsfQbWmxe5GTAvoZ6pOKqChdTlfUXqSNRaeC8YgSOCDWeLRTNfMla2PgPr\nzYvcDJiXUM9UHFXBQuriviJ2JGosvBeMQBHBhpPFopkvGStbn4H15kW0wBUDMLjhmYrjQHAlJxas\nOhP3FbFzSHWZsHMvGDg2nCwiN0ZT2m6hgRDJITDQeSV4carjWTMAkxqenS7Pwue6qDh2B9dycnIF\nU7MN0tZXzJcsvqU2J92XYgUJgZHAguEFSJrbAAACRklEQVQYWQdcMwCTGZ69fDmTomBfxGQ2ZiqO\n3cG1nJZn+G+SFvcVsSPRhhR+914wAkUEG04WiwaRaiYH4YUAQWCw3DiBEky9YFUDMJnhGRmfrfnU\ntlVa9gZXc1pMzTZJi/sKxN9K6mpU05L2ITACFg7J42HRcHAQXggQBAZdBHi9GVnDs24A5tmSNcjr\n/Ky5F2fbNU3ftPNlMWalkclszDQcey+tek7B1MzTq5KW9xWpI9GGEvPuvWAEjgg2nC0WzXzJWNn6\nDKw3L6IAXj4KjBqezQ/nZiKNmYKjM6SZv219eXNmwaueE5ua1UnL+4rUkahFdS8YgSOCDWeLRTNf\nMla2PgPrzYtogKsGYELDs7fTP/r8wmQ2ZhqOvZdWNScyNdsiLe8rEKMsoavRhpphNwIDjeTxcLxm\ndjA4CBACBIERCocVKMHUAtcMwOSGZ+L7tsqLra8QtZwSU7M66Y6+InQkWlcp294LRiCFYMPpYdF6\nnaeYx3oG4YUAQWAkyYHhGFkDXBqAyQ3PvFP6K9NozjQcey+tMqfE1GyDdEdfaeZpB00BU2CtQIfh\nGRufrUF2tp2amm2Ttr6ys7IZncdQIDMAa6VEhmc0ts6987Esp9nUbJO09ZU718qWf0gFMgOwZoZk\neEZj8+S7HsxyCqZmm6Str9y1VLb4gyrwawZgv6hfT07WV36xMLaUKfBHFLC+8kcKbWmaAr+oQOgr\ns5f18y8ua0uZAqbAYyrwPncT988V3OMh/Ov8mGlaVqaAKfCLCnyFdjL9H+4hgujFoAUMAAAAAElF\nTkSuQmCC\n",
      "text/latex": [
       "$$\\left[\\begin{array}{ccccccccccccccc}ig_{0 1} & 0 & 0 & igxp_{11} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\\\0 & ig_{0 2} & 0 & 0 & igxp_{12} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\\\0 & 0 & ig_{0 3} & 0 & 0 & igxp_{13} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\\\xm_{01} & 0 & 0 & c_{11} & zp_{12} & zm_{13} & xp_{21} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\\\0 & xm_{02} & 0 & zm_{11} & c_{12} & zp_{13} & 0 & xp_{22} & 0 & 0 & 0 & 0 & 0 & 0 & 0\\\\0 & 0 & xm_{03} & zp_{11} & zm_{12} & c_{13} & 0 & 0 & xp_{23} & 0 & 0 & 0 & 0 & 0 & 0\\\\0 & 0 & 0 & xm_{11} & 0 & 0 & c_{21} & zp_{22} & zm_{23} & xp_{31} & 0 & 0 & 0 & 0 & 0\\\\0 & 0 & 0 & 0 & xm_{12} & 0 & zm_{21} & c_{22} & zp_{23} & 0 & xp_{32} & 0 & 0 & 0 & 0\\\\0 & 0 & 0 & 0 & 0 & xm_{13} & zp_{21} & zm_{22} & c_{23} & 0 & 0 & xp_{33} & 0 & 0 & 0\\\\0 & 0 & 0 & 0 & 0 & 0 & xm_{21} & 0 & 0 & c_{31} & zp_{32} & zm_{33} & xp_{41} & 0 & 0\\\\0 & 0 & 0 & 0 & 0 & 0 & 0 & xm_{22} & 0 & zm_{31} & c_{32} & zp_{33} & 0 & xp_{42} & 0\\\\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & xm_{23} & zp_{31} & zm_{32} & c_{33} & 0 & 0 & xp_{43}\\\\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & ogxm_{31} & 0 & 0 & og_{41} & 0 & 0\\\\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & ogxm_{32} & 0 & 0 & og_{42} & 0\\\\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & ogxm_{33} & 0 & 0 & og_{43}\\end{array}\\right]$$"
      ],
      "text/plain": [
       "⎡ig₀ ₁    0      0    igxp₁₁    0       0      0     0     0      0       0   \n",
       "⎢                                                                             \n",
       "⎢  0    ig₀ ₂    0      0     igxp₁₂    0      0     0     0      0       0   \n",
       "⎢                                                                             \n",
       "⎢  0      0    ig₀ ₃    0       0     igxp₁₃   0     0     0      0       0   \n",
       "⎢                                                                             \n",
       "⎢xm₀₁     0      0     c₁₁     zp₁₂    zm₁₃   xp₂₁   0     0      0       0   \n",
       "⎢                                                                             \n",
       "⎢  0    xm₀₂     0     zm₁₁    c₁₂     zp₁₃    0    xp₂₂   0      0       0   \n",
       "⎢                                                                             \n",
       "⎢  0      0    xm₀₃    zp₁₁    zm₁₂    c₁₃     0     0    xp₂₃    0       0   \n",
       "⎢                                                                             \n",
       "⎢  0      0      0     xm₁₁     0       0     c₂₁   zp₂₂  zm₂₃   xp₃₁     0   \n",
       "⎢                                                                             \n",
       "⎢  0      0      0      0      xm₁₂     0     zm₂₁  c₂₂   zp₂₃    0      xp₃₂ \n",
       "⎢                                                                             \n",
       "⎢  0      0      0      0       0      xm₁₃   zp₂₁  zm₂₂  c₂₃     0       0   \n",
       "⎢                                                                             \n",
       "⎢  0      0      0      0       0       0     xm₂₁   0     0     c₃₁     zp₃₂ \n",
       "⎢                                                                             \n",
       "⎢  0      0      0      0       0       0      0    xm₂₂   0     zm₃₁    c₃₂  \n",
       "⎢                                                                             \n",
       "⎢  0      0      0      0       0       0      0     0    xm₂₃   zp₃₁    zm₃₂ \n",
       "⎢                                                                             \n",
       "⎢  0      0      0      0       0       0      0     0     0    ogxm₃₁    0   \n",
       "⎢                                                                             \n",
       "⎢  0      0      0      0       0       0      0     0     0      0     ogxm₃₂\n",
       "⎢                                                                             \n",
       "⎣  0      0      0      0       0       0      0     0     0      0       0   \n",
       "\n",
       "    0      0     0     0  ⎤\n",
       "                          ⎥\n",
       "    0      0     0     0  ⎥\n",
       "                          ⎥\n",
       "    0      0     0     0  ⎥\n",
       "                          ⎥\n",
       "    0      0     0     0  ⎥\n",
       "                          ⎥\n",
       "    0      0     0     0  ⎥\n",
       "                          ⎥\n",
       "    0      0     0     0  ⎥\n",
       "                          ⎥\n",
       "    0      0     0     0  ⎥\n",
       "                          ⎥\n",
       "    0      0     0     0  ⎥\n",
       "                          ⎥\n",
       "   xp₃₃    0     0     0  ⎥\n",
       "                          ⎥\n",
       "   zm₃₃   xp₄₁   0     0  ⎥\n",
       "                          ⎥\n",
       "   zp₃₃    0    xp₄₂   0  ⎥\n",
       "                          ⎥\n",
       "   c₃₃     0     0    xp₄₃⎥\n",
       "                          ⎥\n",
       "    0     og₄₁   0     0  ⎥\n",
       "                          ⎥\n",
       "    0      0    og₄₂   0  ⎥\n",
       "                          ⎥\n",
       "  ogxm₃₃   0     0    og₄₃⎦"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "display(InvM)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}
